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1. Scope of tlie Pragroin 



The program is a general purpose finite-difTerence FORTRAN-77 
code designed to solve steady-state three-dimensional heat transfer and 
fluid flow problems in rectangular coordinates. The code is designed to 
effectively handle conjugate domain problems, arising when both solid and 
fluid are present in the computational region. Pure heat conduction 
problems or only fluid flow problems can also be solved. 

A control-volume formulation is used to discretize the governing 
equations expressed in a canonical form. The control volumes for the 
velocity components are staggered with respect to those for the temperature 
and pressure. Harmonic mean formulation for the interface diffusivities is 
used to effectively handle sharp discontinuities in property values in the 
computation domain, without having to resort to very fine mesh in the 
region of discontinuities. The velocity-pressure coupling is handled by the 
SIMPLER algorithm of Patankar [1]. The solution to the steady state 
problem is achieved by providing an initial guess of all the dependent 
variables and proceeding through an iterative scheme using the Thomas 
algorithm (also known as the Tri-Diagonal Matrix Algorithm). Iterations 
are to be stopped when values of variables in successive iterations do not 
change more than a specified convergence criterion. 

2. Program Stnjcture 

The code is highly modular consisting of a short main program and 
various subroutines. The problem independent portion consists of the main 
program and subroutines PROFIL, TDMA, FORM and OPTION. The 
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j)robleni dependent part, wliich the user is primarily concerned with, 
consists of a BLOCK DATA subprogram and a subroutine USE. 

A brief description of the various subprograms/subroutines follows; 



PROFIL 


The profile assumptions in between node-points for the 
dependent variables are incorporated here. Presently, it 
contains the power law profiles [1]. Other profile schemes 
such as the hybrid scheme or the central differencing 
scheme may be incorporated by the user. 


TDMA 


The line-by-line TDMA (Tri-Diagonal Matrix Algorithm) is 
used for solving the equations. The subroutine solves for a 
canonical equation described later (see Eq. 1). 


FORM 


Evaluates the geometric parameters such as lengths, areas 
etc. in ENTRY GEOMET and the coefficients for the various 
equations in ENTRY COEFF which are then used to solve 
for various equations in accordance with the SIMPLER 
algorithm [1]. 


OPTION 


The user may access two optional utilities (entries), 
UMESH and PRINT. UMESH is to be accessed if a uniform 
grid is required, while PRINT is called if dependent 
variables are to be printed out over the entire domain. 
Calling PRINT causes a plane by plane field printout, for 



different planes perpendicular to the z axis. 

BLOCK DATA Contains block COMMONS that contain default values of 
various quantities described subsequently. The default 
values can be changed by the user. 
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USE 



Contains problem dependent information to be specified by 
the user. The different entries in the subroutine USE are 
MESH, BEGIN, VARRIIO, BNDRY, PRTOU’T and DIFFUS. 
The usage of these entries is described later. 

3. Program Usage 



It is required at first to cast the equations to be solved in the following 
canonical form; 



^{pU(})) + -^(pV(t)) + ;^(pW(l)) 



— fr— 1 + —(r—) + — fr— -1 f s 



( 1 ) 



where p is the fluid density, u, v and w are the velocity components in the k , 
y and z directions respectively, ({) is the dependent variable to be solved for, T 
is the appropriate diffusivity and S is the source term. 

For efficient execution of the algorithm, and for better convergence 
pi'ospects, it is necessary that S be represented in the form 

S=Sq+Sp<|) (2) 

and it is required that Sp < 0. For an appropnate linearization of the 
source term S, when it does not have a readily obtainable form of Eq. (2), the 
user is referred to Patankar [1], pp. 48-49. 
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It is not required to cast the continuity equation into the canonical 
form, since this is incorporated in the program internally. For example, 
when solving a combined incompressible fluid-flow, heat transfer problem, 
the user needs to cast the three momentum equations and one energy 
equation in the canonical form of Eq. (1). Note that the pressure gradient 
terms encountered in the momentum equations are incorporated internally 
and the user is cautioned not to incorporate them in the source term S. 

As mentioned earlier, the user needs to be concerned only with the 
BLOCK DATA and subroutine USE. A detailed description of these 
segments now follows. 

BLOCK DATA 

The array F(I,J,K,NF) contains different variables or properties. The 
default equivalence of the array F is given below: 



NF 



Quantity Stored in F(l.J.K.NF) FORTRAN Array in USE 
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2 



3 



1 



velocity in x direction, u 
velocity in y direction, v 
velocity in z direction, w 
pressure correction 



V(I,J,K) 

W(I,J,K) 



U(I,J,K) 



PC(I,J,K) 
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6 



pressure 



temperature 



T(I,J,K) 

P(I,J,K) 
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The BLOCK DATA contains several data that nuay bo reinitialized by 
the user. Each of these is explained below; 



IPREF, JPREF, KPREF 


Values of the indices (I,J,K) where the pressure is 
set to zero. Note that the program evaluates only 
pressure differences and not the actual pressures. 


NTIMES - 


The value of NTIMES(NF) determines the number 
of sweeps of the TDMA solution over the entire 
domain during one iteration. More than one 
sweep per iteration may be required for higlily 
nonlinear equations to assure convci'gcnce. 
Clearly, a larger value of NTIMES(NF) will result 
in a longer computational time. 


LAST - 


Defines the total number of iterations to be 
performed. 


LPRINT - 


The default value of this logical FORTRAN 
variable is FALSE. Setting LPRINT(NF) to be true 
will cause a field printout of the particular 
variable associated with NF when a call to the 
entry PRINT is invoked through the statement 
CALL PRINT. 


LSOLVE - 


The default value of this logical FORTRAN 
variable is FALSE. Setting LSOLVE(NF) to be true 
causes the solution of the particular variable 
associated with NF. In a pure conduction 
problem, for example, only LSOLVE(5) must be set 
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RELAX - 



RHOCON - 
TITLE(NF) - 



true. In a pure fluid flow problem (no lieat 
transfer) LSOLVE(NE), NF=1,2,3,4, and G must ho 
declared true. In a heat transfer and fluid flow 
problem, LSOLVE(NF), NF=1 to 6 should bo set 
true. 

The value for RELAX(NF) is the underrelaxation 
factor for the solution of the variable associated 
with NF. RELAX(NF)>1 causes overrelaxation, 
while RELAX(NF)<1 causes underrelaxation. 
Underrelaxation may be necessary in many 
problems when convergence is difficult to obtain, 
and will increase the computational time. 

Value of p in Eq.(l) is specified here. 

This describes the title to be given to each field 
printout. 



Subroutine USE contains various segments (entries) called by the 
main program. The purpose of these and their usage is described here. 

See Appendix A for meaning of the various symbols. 

MESH: The user supplies a mesh in this entry. The computational 

domain must be divided into rectangular control volumes. The 
internal grid points are placed at the geometric centers of the 
control volumes. The boundary grid points are placed on the 
geometric centers on the control volumes faces that form part 
of the boundary (See Figs. A-1 and A-2 in Appendix A). 
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BEGIN: 



Two options are available- 

1. Uniform grid - The user supplies values of LI, Ml and Nl; 
XL, YL and ZL (see Appendix A); and then invokes CALL 
UMESH to form a uniform mesh, i.e. all the control volumes in 
a particular direction have the same length; 

ENTRY MESH 

Ll=10 

Ml=10 

Nl=20 

XL=2. 

YL=3. 

ZL=1. 

CALL UMESH 
RETURN 

This implies a computational domain that is 2, 3 and 1 units 
long in the x, y and z directions respectively. All the control 
volumes thus have dimensions of 2/10, 3/10 and 1/20 units in 
the X, y and z directions respectively. 

2. Non-uniform grid - The user supplies values of (XU(I), 
I=1,L1), (W(J), J=1,M1), ZW(K), K=1,N1); XL, YL, ZL; LI, Ml 
and Nl. The code internally determines the locations of the 
boundary and internal grid points (stored in X(I),Y(J),Z(K)) 
which are placed at the geometric centers of the control 
volumes. For example, X(I)=0.5*[XU(IHXU(N1)] etc. 
Boundary points are grid points that are located on planes x=0 
or x=XL or y=0 or y=YL or z=0 or z=ZL. 

The user specifies/calculates quantities that are required prior 
to the start of iterations such as parameter values, initial 
guess to the solution and boundary conditions. 
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Note: 



Entries MESTI and liEGIN are called by the main program 



VARRHO: 

BNDRY: 

Note: 



only once before the iterations begin. 



Functional dependency of the fluid density can be accounted for 
in this entry for variable density situations. 



The boundary values of the variables are updated every 
iteration in accordance with the boundary conditions. 

Boundary conditions that involve specification of the dependent 
variable value can be incorporated in BEGIN and it is not 
necessary to include these in BNDRY. This is so because the 
program solves for the dependent variables only at the internal 
grid points and not at the boundary points. 

Examples of different boundary conditions: 



ENTRY BNDRY 

DO 1 1=1, LI 
DO 1J=1,M1 

C Insulated boundary 

C At z=ZL, dT/dz=0 

T(I,J,N1)=T(I,J,N1-1) 

C Convective condition 

C At z=0, h(Tf-T)=-kdT/dz 

C H=h (convection coefficient), TF=Tf (fluid temperature) 

C TK=k (thermal conductivity) 

DELZ=Z(2)-Z(1) 

T(I,J,1)=(H*DELZ*TF+TK*T(I,J,2))/(TK+H''DELZ) 

C Specified heat flux condition 

C at z=ZL, -kdT/dz=-q", q">0 

C QPRIME=q" (heat flux) 

DELZ=Z(N1)-Z(N1-1) 

T(I,J,N1)=T(I,J,N1-1)+QPRIME*DELZ/TK 
1 CONTINUE 
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RETURN 



PRTOUT: 



DIFFUS: 



It may be necessary to monitor values of different variables or 
other quantities of interest as each iteration proceeds. This 
can be done in ENTRY PRTOUT. At the end of the iterations, 
CALL PRINT invokes a field printout of the variables for which 
LPRINT has been set true. 



The user specifies P, and Sp (Eqs. 1,2) in the variables 

CAM(I,J,K), CON(I,J,K) and AP(I,J,K), respectively. The 
diffusivity P is specified at all points, while and Sp are 

specified at the internal points (i.e. all points excluding the 
boundary points). For example: 



ENTRY DIFFUS 

C Specify NF since DIFFUS is called by the main program 
C more than once every iteration. 

DO 2 I=1,L1 
DO 2 J=1,M1 
D02K=1,N1 

C Vise is the viscosity 

IF((NF.EQ.l).OR.(NF.EQ.2).OR.(NF.EQ.3))GAM(I,J,K)= 

IVISC 

C TK is the thermal conductivity, CP is the specific heat 
C capacity. 

IF(NF.EQ.5)GAM(I,J,K)=TK/CP 
2 CONTINUE 

ENDIF 

C In the event of uniform, volumetric heat generation in 
C the domain (Q W/m^). 

DO 3 I=2,L1-1 
D0 3 J=2,M1-1 
DO 3 K=2,N1-1 

IF(NF.EQ.5)CON(I,J,K)=Q/CP 
C If AP(I,J,K) is not specified, it is zero by FORTRAN 
C default. 
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CONTINUE 



C SPECIAL TREATMENT OF FLUX TYPE BOUNDARY 
C CONDITIONS CrO IMPROVE CONVERGENCE) 

C SET BOUNDARY VAI.UE OF DIFFUSIVITY TO ZERO. 
C CONVERT SURFACE HEAT FLUX TO AN 
C EQUIVALENT HEAT SOURCE IN THE CONTROL 
C VOLUME ADJACENT TO THE BOUNDARY. 

DO 4 1=1, LI 
D0 4 J=1,M1 

C Insulated boundary 

C At z=ZL, dT/dz=0, 

IF(NF.EQ.5)GAM(I,J,N1)=0.0 
C Convective condition 

C At z=0, h(T(-T)=-kdT/dz 

C S=Equivalent heat source/volume=surface heat Ilu.x '' 

C surface area/volume (source term for next to boundary 
C node). 

C Surface convective flux=h(Tf-T)=(Tf-Ti)/[l/li4-(zi-z|^)/k] 

C where Tj is the temperature at the node adjacent 
C to the boundary, z[ is the z location of the node, zjj is the 

C z location of the boundary. The factor in square brackets 

C is the thermal resistance (Rth) to convection and 

C conduction. 

C Surface area=XCV(I)*YCV(J). 

C Volume of control volume=XCV(I)*YCV(J)*ZCV(2). 

C Surface area/volume=l./ZCV(2) 

C Ti = T(I,J,2), zi-zb=Z(2)-Z(l), T=T(I,J,1). 

C Source S may be split into S^ and Sp 
C Sc=(Tf/Rfh)/ZCV(2), SD=-l/(Rth*ZCV(2)) 

RTH=(l.ffl+(Z(2)-Z(l))^K) 

IF(NF.EQ.5)THEN 

GAM(I,J,1)=0.0 

C Add Sc and Sp for the control volume adjacent to the 

C boundary. 

CON(I,J,2)=CON(I,J,2)+TF/(RTH*ZCV(2)) 

AP(I,J,2)=AP(I,J,2)-1./(RTH*ZCV(2)) 

ENDIF 

C Specified heat flux condition 

C At z=ZL, -kdT/dz=-q", q">0 

C QPRIME=q" (heat flux) 

IF(NF.EQ.5)THEN 

GAM(I,J,N1)=0.0 

CON(I,J,Nl-l)=CON(I,J,Nl-l)+QPRIME/ZCV(Nl-l) 

ENDIF 

4 CONTINUE 

RETURN 
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Conjugate Domains Frequently it is required to calculate heat transfer 
in domains that contain both fluid and solid. The solid then can be easily 
modeled by setting the fluid viscosity to a very large value (eg. lO^O) in the 
solid region of the computational domain. The velocities in that region then 
will be computed by the code to be zero. Appropriate solid properties such 
as thermal conductivity can be prescribed in the solid region. Care must be 
taken that the control volume faces coincide with the actual physical 
discontinuities in the domain. 

4. Program Flow-chart 

The flow-chart for the program execution is shown in Fig. 1. Tlie 
usage and the purpose of the various subroutines and subprograms has 
already been explained in the preceding sections. A single iteration 
consists of execution of the SIMPLER algorithm to evaluate the velocity and 
pressure field in the domain and then the solution of any other dependent 
variables to be solved for, based on the calculated velocities. 
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Fig. 1 Program Flow-chart. 
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y\PPENDIX A - OUier FORTRAN Variables in USE 
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ITER 



L1,M1,N1 

XL,YL,ZL 

X(I),Y(J),Z(K) 

XU(I),YV(J),ZW(K) 



FX(I), FXM(I) 



Counter for iterations. Program stops when 
ITER=LAST. 

Number of grid points in the X, Y and Z directions, 
respectively. Ll-2, Ml-2 and Nl-2 are the total 
number of control volumes in the X, Y and Z 
directions. 

Domain lengths in the X,Y, and Z directions, 
respectively. 

Locations of the node points where temperatures and 
pressures are stored. The temperature T(I,J,K) is 
stored at the location (X(I),Y(J),Z(K)). Note 
X(1)=0.,X(L1)=XL. 

Locations of control volume faces. XU(I) is the 
location of the control volume face perpendicular to 
the X axis for the control volume around the point 
(X(I),Y(J),Z(K)) such that XU(D<X(I). 
XU(2)=0.,XU(L1)=XL. XU(1) is meaningless. 

U(I,J,K) is evaluated at the point (XU(I),Y(J),Z(K)). 
Similar description is applicable for YV(J) and 
ZW(K). 

Interpolation factors in the x direction that enable 
calculation of a quantity a control volume interface 
perpendicular to the x direction. EXAMPLE: 
Temperature at the interface of two control volumes, 
say at XU(I), can be calculated as 
FX(I)*T(I,J,K)+FXM(I)*T(I-1,J,K) 
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XDIF(I),YDIF(J), 

ZDIF(K) 

XCV(I), YCV(J), 
ZCV(K) 



Similar interpretation for FY(J), FYM(J); FZ(K), 
FZM(K). 

XDIF(I)=X(I)-X(I-1), YDIF(J)=Y(J)-Y(J- 1), 
ZDIF(K)=Z(K)-Z(K-1) 



Dimensions in the x, y and z directions respectively 
associated with control volume around the point 
(X(I),Y(J),Z(K)). eg. XCV(I)=XU(I+1)-XU(I). Thus 
XCV(l) and XCV(L1)=0. 
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(XL.YL.O) 



(O.YL.O) 




(XL.YL.ZL) 



Ml 



I control volumes 



Fig. A-1. The computational domain. The coordinate system is right-handed. 
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Plane ZW(K) 




here. 



Fig. A-2. The control volume. 
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APPENDIX B - Program Source Code. The subroutine USE is for the 
problem described in Appendix C. 
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ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c 

C A GENERAL PURPOSE FORTRAN PROGRAM FOR SOLVING THREE- 

C DIMENSIONAL HEAT TRANSFER AND FLUID FLOW IN RECTANGULAR 

C COORDINATES 

C 

C VERSION 1.00, JULY 1990 

C 

C DEVELOPED BY: 

C S. B. SATHE 

C CODE 69ST, department OF MECHANICAL ENGINEERING 

C U.S. NAVAL POSTGRADUATE SCHOOL 

C MONTEREY, CA 93943 

C TEL: ( 408 )-646-2417 

C BITNET ADDRESS: 5140P@NAVPGS 

C 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

(]]**********A******:fc****AA******A********************AA***A 

LOGICAL LSTOP 
COMMON/CNTL/LSTOP 

(^AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 



CALL 


MESH 


CALL 


GEOMET 


CALL 


BEGIN 


10 CALL 


VARRHO 


CALL 


BNDRY 


CALL 


PRTOUT 


IF(LSTOP) STOP 


CALL 


COEFF 



GO TO 10 
END 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

SUBROUTINE PROFIL 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

COMMON/COEF/FLOW , DI FF , ACOF 

^AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 



ACOF=DI FF 

IF( FLOW. EQ. 0 . ) RETURN 
TEMP=DIFF-ABS ( FLOW) *0.1 
ACOF=0 . 

IF(TEMP.LE.O. ) RETURN 
TEMP=TEMP/DIFF 
ACOF=DIFF*TEMP**5 
RETURN 



END 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

SUBROUTINE TDMA 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

LOGICAL LSOLVE, LPRINT, LBLK, LSTOP 

COMMON F(17,17,13,5),P(17,17,13) , RHO ( 17 , 1 7 , 1 3 ) ,GAM(17 
1 CON( 17,17,13),AKP(17,17,13) , AKM( 17, 17, 13), AP (17, 17 

17,13) , AJP( 17 , 17 



AIP( 17 
COMMON 
X( 17 ) , XU 



17,13) , AIM( 17 



13 ) , AJM( 17 



17 
13) , 
17,13) 



17 ) ,XDIF( 17 ) ,XCV 
YDIF( 17 ) ,YCV 



17 ) ,XCVS( 17 ) 
17 ) , YCVS( 17 ) 



Y( 17 ) , YV( 17 ) 

Z(17),ZW(17),ZDIF(17),ZCV(17) , ZCVS( 17 ) , 

YCVR( 17 ) , YCVRS( 17 ) , ARX{ 17 ) , ARXJ( 17 ) , ARXJP( 17 ) , 

R( 17 ) , RMN( 17 ) , SX( 17 ) , SXMN( 17 ) ,XCVI ( 17 ) , XCVIP( 17 ) 
YCVJ ( 17 ) , YCVJP( 17 ) , ZCVK( 17) ,ZCVKP(17) 



13) , 
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COMMON DU (17, 17, 13), DV (17, 17, 13), DW (17, 17, 13), FV (17), FVP( 17 ) , 

1 FX(17) ,FXM(17) ,FY(17) ,FYM(17) , PT ( 1 7 ) , QT ( 1 7 ) , 

2 FZ ( 17 ) , FZM( 17 ) 

COMMON/ I NDX/RELAX( 13),LPRINT(13),LBLK(11),NTIMES(10), 

1LS0LVE( 10 ) , TIME, DT, XL, YL , ZL, S , RHOCON, ZERO, 

2NF,NFMAX,NP, NRHO, NOAM, LI , L2 , L3 , Ml , M2 , M3 , N1 , N2 , N3 , 

3IST, JST, KST, ITER, LAST, 

4ipref,jpref, kp ref, mode 

DIMENSION D( 17 ) , VAR( 17 ) , VARM( 17 ) , VARP( 17 ) , PH I BAR ( 17 ) 

COMMON/HEADIN/TITLE 

CHARACTER*10 TITLE(13) 

ISTF=IST-1 

JSTF=JST-1 

KSTF=KST-1 

IT1=L2+IST 

IT2=L3+IST 

JT1=M2+JST 

JT2=M3+JST 

KT1=n2+KST 

KT2=N3+KST 



NTSSS=NTIMES(NF) 

DO 999 NT=1,NTSSS 
DO 391 N=NF,NF 
C 

IF( .NOT.LBLK(NF) )GO TO 60 
60 CONTINUE 

COMMENCE TDMA LINE BY LINE SWEEPS FOR SOLUTION 
DO 90 K=KST,N2 
DO 90 J=JST,M2 
PT( ISTF)=0 . 

QT(ISTF)=F( ISTF, J,K,N) 

DO 70 I=IST,L2 

DENOM=AP( I,J,K)-PT(I-1)*AIM(I,J,K) 

PT( I )=AIP( I , J , K)/DENOM 

TEMP=CON( I , J, K)+AJP( I , J, K) *F( I , J+ 1 , K , N ) +AJM ( I,J,K)*F(I,J-1,K,N) 
1 +AKP( I,J,K)*F(I,J, K+1 ,N)+AKM( I,J,K)*F(I,J,K-1,N) 

QT( I )=( TEMP+AIM( I,J,K)*QT(I-1) )/DENOM 
70 CONTINUE 

DO 80 II=IST,L2 
I=IT1-II 

80 F(I,J,K,N)=F(I+1,J,K,N)*PT(I )+QT( I ) 

90 CONTINUE 



DO 190 KK=KST,N3 
K=KT2-KK 

DO 190 JJ=JST,M3 

J=JT2-JJ 

PT( ISTF)=0 . 

QT( ISTF )=F( ISTF, J , K,N ) 

DO 170 I=IST,L2 

DENOM=AP (I,J,K)-PT(I-1)*AIM(I,J,K) 

PT( I )=AIP( I , J, K)/DENOM 

TEMP=CON( I , J, K)+AJP( I , J , K ) * F ( I , J+ 1 , K , N ) +AJM ( I , J, K) *F( I , J-1 , K,N) 
1 +AKP( I,J,K)*F(I,J, K+1 ,N )+AKM( I,J,K)*F(I,J,K-1,N) 

QT( I )=(TEMP+AIM( I,J,K)*QT(I-1) )/DENOM 
170 CONTINUE 

DO 180 II=IST,L2 
I=IT1-I I 
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180 F(I,J,K,N)=F(I+1,J,K,N)*PT(I )+QT( I ) 
190 CONTINUE 



DO 290 I=IST,L2 
DO 290 K=KST,N2 
PT( JSTF)=0 . 

QT(JSTF)=F(I,JSTF,K,N) 

DO 270 J=JST,M2 

DENOM=AP( I,J,K)-PT(J-1)*AJM(I,J,K) 

PT ( J ) =AJ P ( I , J , K ) /DENOM 

TEMP=CON( I , J , K) +AKP( I,J,K)*F(I,J,K+1 , N) +AKM( I,J,K)*F(I,J,K-1,N) 
1 +AIP(I,J,K)*F(I+1,J, K ,N)+AIM( I,J,K)*F(I-1,J,K,N) 

QT( J )=(TEMP+AJM( I,J,K)*QT(J-1) ) /DENOM 
270 CONTINUE 

DO 280 JJ=JST,M2 
J=JTl-JJ 

280 F(I,J,K,N)=F(I,J + 1,K,N) *PT( J)+QT( J) 

290 CONTINUE 



DO 390 II=IST,L3 
I=IT2-I I 

DO 390 KK=KST,N3 

K=KT2-KK 

PT( JSTF ) =0 . 

QT( JSTF)=F( I , JSTF,K,N) 

DO 370 J=JST,M2 

DENOM=AP( I,J,K)-PT(J-1) *AJM( I , J , K) 

PT ( J ) =AJ P ( I , J , K ) /DENOM 

TEMP=CON( I , J , K )+AKP( I,J,K)*F(I,J,K+1,N) +AKM( I,J,K)*F(I,J,K-1,N) 
1 +AIP(I,J,K)*F(I+1,J,K,N)+AIM(I,J,K)*F(I-1,J,K,N) 

QT( J )=( TEMP+AJM( I,J,K)*QT(J-1) ) /DENOM 
370 CONTINUE 

DO 380 JJ=JST,M2 
J=JTl-JJ 

380 F{I,J,K,N)=F(I , J+1 , K,N) *PT( J )+QT( J ) 

390 CONTINUE 



DO 490 J=JST,M2 
DO 490 I=IST,l2 
PT( KSTF)=0 . 

QT(KSTF)=F(I,J,KSTF,N) 

DO 470 K=KST,N2 

DENOM=AP( I , J , K )-PT( K-1 ) *AKM( I , J , K) 

PT( K )=AKP( I , J , K ) /DENOM 

TEMP=CON( I,J,K)+AIP(I,J,K)*F(I+1,J,K,N)+AIM(I,J,K)*F(I-1,J,K,N) 
1 +AJP( I,J,K)*F(I,J+1,K,N) +AJM( I,J,K)*F(I,J-1,K,N) 

QT( K )=( TEMP+AKM( I , J , K ) *QT( K-1 ) ) /DENOM 
470 CONTINUE 

DO 480 KK=KST,N2 
K=KTl-KK 

480 F(I,J,K,N)=F(I,J,K+1,N) *PT{ K)+QT( K) 

490 CONTINUE 



DO 590 JJ=JST,M3 
J=JT2-JJ 

DO 590 II=IST,L3 

I=IT2-II 

PT( KSTF)=0 . 

QT(KSTF)=F( I ,J,KSTF,N) 
DO 570 K=KST,N2 
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DENOM=AP( I , J,K )-PT( K-1 ) *AKM( I , J, K) 

PT(K)=AKP( I, J,K)/DENOM 

TEMP=CON( I , J , K) +AIP( I,J,K)*F(I+1,J,K,N) +AIM( I,J,K)*F(I-1,J,K,N) 
1 +AJP( I,J,K)*F(I,J+1, K,N) +AJM( I,J,K)*F(I,J-1,K,N) 

QT( K )=( TEMP+AKM( I , J , K ) *QT( K-1 ) )/DENOM 
570 CONTINUE 

DO 580 KK=KST,N2 
K=KT1-KK 

580 F(I,J,K,N)=F(I,J,K+1,N)*PT( K) +QT( K ) 

590 CONTINUE 

C 

391 CONTINUE 

^ic-k'kiKif-k-k-k'k-kif-k-k-k-kiK-k-k^if-k-k-kiK-k-k'kii'k-k-k'k'M-k^-k-k'kit'k'k'k'k'kif-k'k'k'k'k-k'k'k'k'k-k'k 

999 CONTINUE 

*********************************************** 

ENTRY RESET 

Q*********************************************** 

DO 400 K=2,n2 
DO 400 J=2,M2 
DO 400 1=2, L2 
CON( I, J,K)=0. 

AP( I, J,K)=0. 

400 CONTINUE 
RETURN 
END 
C 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

SUBROUTINE FORM 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

LOGICAL LSOLVE,LPRINT,LBLK,LSTOP 

COMMON F(17,17,13,5),P(17,17,13) ,RHO( 17,17,13) ,GAM( 17, 17 , 13 ) , 

1 CON( 17 , 17,13) ,AKP( 17,17,13) ,AKM( 17, 17, 13), AP (17, 17, 13), 

2 AIP(17,17,13),AIM(17,17,13),AJP(17,17,13) , AJM( 17 , 17 , 13 ) 

COMMON 

3 X(17),XU(17),XDIF(17) ,XCV( 17 ) ,XCVS( 17 ) , 

4 Y(17),YV(17),YDIF(17), YCV( 17 ) , YCVS ( 17 ) , 

5 2(17), ZW (17), 2DIF( 17 ) , ZCV( 17 ) , 2CVS( 17 ) , 

6 YCVR( 17 ) , YCVRS ( 17 ) , ARX( 17 ) ,ARXJ ( 17 ) , ARXJP( 17 ) , 

7 R( 17 ) , RMN( 17 ) , SX( 17 ) , SXMN( 17 ) ,XCVI ( 17 ) , XCVIP( 17 ) , 

8 YCVJ ( 17 ) , YCVJP ( 17 ) , ZCVK( 17 ) , ZCVKP( 17 ) 

COMMON DU (17, 17, 13), DV (17, 17, 13), DW (17, 17, 13), FV (17), FVP (17), 

1 FX( 17 ) , FXM( 17 ) , FY( 17 ) , FYM( 1 7 ) , PT ( 1 7 ) , QT ( 1 7 ) , 

2 FZ( 17) ,FZM( 17) 

COMMON/I NDX/RELAX( 13 ) , LPRINT( 13 ) ,LBLK( 11 ) ,NTIMES( 10 ) , 

1LS0LVE( 10) ,TIME,DT,XL,YL,ZL,S,RHOCON,ZERO, 

2NF,NFMAX,NP, NRHO,NGAM, LI , L2 , L3 , Ml , M2 , M3 , Nl , N2 , N3 , 

3 1 ST, J ST, KST, ITER, LAST, 

4IPREF, JPREF, KPREF,MODE 
COMMON/HEAD I N/TI TLE 
CHARACTER*10 TITLE(13) 

COMMON/CNTL/LSTOP 
COMMON/SORC/SMAX , SSUM 
COMMON/COEF/FLOW , DI FF , ACOF 

DIMENSION U(17,17,13),V(17,17,13),W(17,17,13),PC(17,17,13) 
DIMENSION T( 17 , 17 , 13 ) 

EQUIVALENCE( F( 1, 1 ,1,1),U(1,1,1)),(F(1,1,1,2),V(1,1,1)), 

1 (F(1,1,1,3),W(1,1,1)),(F(1,1,1,4),PC(1,1,D) 

2 , (F(l, 1,1,5) ,T(1, 1,1) ) 

DIMENSION C0F1(17,17,13,7) , COF2 ( 17 , 17 , 1 3 , 7 ) , COF 3 ( 1 7 , 1 7 , 1 3 , 7 ) 
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DIMENSION COF4(17,17,13,7) , CONI ( 1 7 , 1 7 , 1 3 ) , CON2 ( 1 7 , 1 7 , 1 3 ) 
DIMENSION CON3 ( 17 , 17 , 13 ) 

1 FORMAT( 15X, 'COMPUTATION IN CARTESIAN COORDINATES') 

2 format! 15X, ' COMPUTATION FOR AXISYMMETRIC SITUATION') 

3 FORMAT( 15X, ' COMPUTATION IN POLAR COORDINATES') 

4 FORMAT( 14X, 40( IH* ) ,//) 

OME HERE TO CALCULATE GRIDS SPECIFICATION 
ENTRY GEOMET 

■k-k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k-k-k-k-k-^-k-k-k-k-k-k-k-k-k-k-kiK-k-k-k-k'k-k-k'k-k'k'k'kit'k'k-k'k'k'k'k 



L2=L1-1 
L3=L2-1 
M2=M1-1 
m3=M2-1 
N2=N1-1 
N3=N2-1 
X( 1 )=XU( 2 ) 

DO 5 1=2, L2 

5 X( I ) =0 . 5* { XU( I+l )+XU( I ) ) 
X( Ll ) =XU( Ll ) 

Y( 1 )=YV( 2 ) 

DO 10 J=2,M2 

10 Y( J )=0 . 5* ( YV( J + 1 )+YV( J ) ) 
Y(M1 )=YV(Ml ) 

Z ( 1 ) =ZW( 2 ) 

DO 7 K=2,N2 

7 Z ( K ) = 0 . 5* ( ZW( K+1 ) +ZW( K ) ) 
Z(N1)=ZW(N1) 



DO 15 1=2, Ll 
15 XDIF( I )=X( I )-X( I-l ) 

DO 18 1=2, L2 
18 XCV( I )=XU( I+l )-XU{ I ) 

DO 20 1=3, L2 
20 XCVS( I )=XDIF{ I ) 

XCVS ( 3 )=XCVS{ 3 )+XDI F( 2 ) 
XCVS(L2)=XCVS(L2)+XDIF(L1) 
DO 22 1=3, L3 
XCVI ( I )=0 . 5*XCV( I ) 

22 XCVIP( I )=XCVI ( I ) 

XCVIP( 2)=XCV( 2) 

XCVI ( L2 ) =XCV( L2 ) 

DO 175 K=2,Nl 
175 ZDIF( K)=Z( K)-Z( K-1 ) 

DO 178 K=2,N2 
178 ZCV( K)=ZW( K+1 )-ZW( K) 

DO 270 K=3,N2 
270 ZCVS( K)=ZDIF( K) 

ZCVS( 3 )=ZCVS( 3 )+ZDIF{ 2 ) 
ZCVS(N2 )=ZCVS(N2)+ZDIF(N1 ) 
DO 272 K=3,N3 
ZCVK( K )=0 . 5*ZCV( K ) 

272 ZCVKP(K)=ZCVK(K) 

ZCVKP( 2)=ZCV( 2) 
ZCVK(N2)=ZCV(N2) 

DO 35 J=2,M1 
35 YDI F( J )=Y( J )-Y( J-1 ) 

DO 40 J=2,M2 
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40 YCV( J)=YV( J+1)-YV( J) 

DO 45 J=3,M2 
45 YCVS( J)=YDIF( J) 

YCVS( 3 ) =YCVS( 3 ) +YDIF( 2 ) 
YCVS(M2)=YCVS(M2)+YDIF(M1 ) 

DO 277 J=3,M3 
YCVJ( J)=0. 5*YCV( J) 

277 YCVJP( J)=YCVJ( J) 

YCVJP( 2 )=YCV( 2 ) 

YCVJ(M2)=YCV(M2) 

IF(MODE.NE. 1 ) GO TO 55 
DO 52 

RMN( J)=l .0 
52 R(J)=1.0 
GO TO 56 

55 DO 50 J=2,M1 

50 R( J)=R( J-1 )+YDIF( J) 

RMN( 2 )=R( 1 ) 

DO 60 J=3,M2 

60 RMN( J ) =RMN( J-1 ) +YCV( J-1 ) 

RMN( Ml ) =R( Ml ) 

56 CONTINUE 

DO 57 J=l,Ml 
SX( J ) =1 . 

SXMN( J)=l . 

IF(MODE.NE. 3) GO TO 57 
SX( J ) =R( J ) 

IF(J.NE.l) SXMN( J)=RMN( J) 

57 CONTINUE 

DO 62 J=2,M2 
YCVR( J ) =R( J ) *YCV( J ) 

ARX( J ) =YCVR( J ) 

IF(MODE.NE.3) GO TO 62 
ARX ( J ) =YCV( J ) 

62 CONTINUE 

DO 64 J=4,M3 

64 YCVRS (J)=0.5*(R(J)+R(J-1))*YDIF(J) 
YCVRS( 3 )=0 . 5* ( R( 3 )+R( 1 ) ) *YCVS ( 3 ) 
YCVRS (M2)=0.5*(R(M1)+R(M3) ) * YCVS ( M2 ) 
IF(MODE.NE.2) GO TO 67 

DO 65 J=3,M3 

ARXJ( J )=0 . 25* ( 1 .+RMN( J )/R( J ) ) *ARX( J ) 

65 ARXJP( J)=ARX( J)-ARXJ( J) 

GO TO 68 

6 7 DO 66 J= 3 , M3 
ARXJ( J)=0.5*ARX( J) 

66 ARXJP( J)=ARXJ( J) 

68 ARXJP( 2 ) =ARX( 2 ) 

ARXJ(M2)=ARX(M2) 

DO 70 J=3,M3 

FV( J )=ARXJP( J ) /ARX( J ) 

70 FVP( J)=l .-FV( J) 

DO 85 1=3, L2 

FX( I )=0 . 5*XCV( I-l )/XDIF( I ) 

85 FXM( I )=1 .-FX( I ) 

FX( 2 )=0 . 

FXM( 2 ) =1 . 

FX( LI )=1 . 

FXM( LI )=0. 

DO 90 J=3,M2 
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FY( J )=0 . 5*YCV( J-1 )/YDI F( J ) 

90 FYM( J )=1 .-FY( J ) 

FY( 2 )=0 . 

FYM( 2 )=1 . 

FY( Ml ) =1 . 

FYM( Ml ) =0 . 

DO 87 K=3,N2 

FZ(K)=0.5*ZCV(K-1)/ZDIF(K) 

87 FZM( K ) =1 . -FZ ( K ) 

FZ( 2 )=0 . 

FZM( 2 ) = 1 . 

FZ(Nl)=l. 

FZM( N1 )=0 . 

CON, AP , U, V, RHO , PC AND P ARRAYS ARE INITIALIZED HERE 
DO 95 K=1,N1 
DO 95 J=l,Ml 
DO 95 1=1, LI 
PC( I,J,K)=0. 

U( I , J , K )=0 . 

V(I,J,K)=0. 

W( I , J , K )=0 . 

CON( I , J , K)=0 . 

AP( I , J , K)=0 . 

RHO( I , J, K)=RHOCON 
P( I, J,K)=0. 

95 CONTINUE 

IF(MODE.EQ.l) PRINT 1 
IF( MODE . EQ. 2 ) PRINT 2 
IF( MODE . EQ. 3 ) PRINT 3 
PRINT 4 
RETURN 
C 

COME HERE TO CALCUALTE COEFFICIENTS FOR FINITE DIFF. EQNS. 

ENTRY COEFF 

QAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

c 

COEFFICIENTS FOR THE U EQUATION 

c 

CALL RESET 
NF=1 

IF( .NOT.LSOLVE(NF) ) GO TO 100 

IST=3 

JST = 2 

KST=2 

CALL DIFFUS 
REL=1 . -RELAX( NF ) 

DO 103 1=3, L2 
DO 103 J=2,M2 
DO 103 K=2,N2 

COEFFICIENTS AEAST AND AWEST 

FL=U( I,J,K)*(FX(I) *RHO( I , J, K ) + FXM( I ) *RHO( I-l , J , K ) ) 

FLP=U( I + 1,J,K)*(FX(I + 1) *RHO( I + l , J , K ) + FXM( I + l ) *RHO( I , J , K ) ) 
FLOW=YCV( J )*ZCV(K)*0.5*(FL+FLP) 

DIFF=YCV( J ) *ZCV( K ) *GAM( I , J , K )/XCV( I ) 

CALL PROFIL 

AIM(I+1,J,K) =ACOF+AMAXl (ZERO, FLOW) 

AIP(I,J,K)=AIM( I+l , J , K ) -FLOW 
COEFFICIENTS ANORTH AND ASOUTH 

FL=XCVI (I)*V(I,J + 1,K)*(FY(J + 1) *RHO( I , J + 1 , K)+FYM( J + 1 ) *RHO( I , J , K ) ) 
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FLM=XCVIP( I-1)*V(I-1,J+1,K)*(FY(J+1) *RHO( I-1,J+1,K)+FYM(J+1)* 

1 RHO( I-l , J, K) ) 

GM=GAM( I , J , K ) *GAM( I , J+1 , K ) 

1 /( YCV( J ) *GAM( I , J + 1 , K ) +YCV( J + 1 ) *GAM( I , J , K ) + 

2 1 . OE-30 ) *XCVI ( I ) 

GMM=GAM( I-l , J , K ) *GAM( I-l , J + 1 , K) 

1 /( YCV( J ) *GAM( I-l , J+1 , K )+YCV( J+1 ) * 

2 GAM( I-l,J,K)+l.E-30) *XCVIP( I-l ) 

DIFF=ZCV(K)*2.*( GM+GMM) 

FLOW=ZCV( K ) * ( FL+FLM) 

CALL PROFIL 

AJM( I , J+1 , K )=ACOF+AMAXl ( ZERO, FLOW) 

AJP( I , J , K )=AJM( I , J+1 , K)-FLOW 
COEFFICIENTS AIN AND AOUT 

FL=XCVI ( I ) *W( I , J , K+1 ) * ( FZ( K+1 ) *RHO( I,J,K+l)+FZM(K+l)*RHO(I,J,K) ) 
FLM=XCVIP( I-l ) *W( I-l , J,K+1 ) *( FZ( K+1 ) *RHO( I-l ,J, K+1 )+FZM( K+1 ) * 

1 RHO( I-l , J , K ) ) 

GM=GAM( I , J , K ) *GAM( I , J , K+1) 

1 /( ZCV( K) *GAM( I , J , K+1 )+ZCV( K+1 ) *GAM( I , J , K) + 

2 1 . OE-30 ) *XCVI ( I ) 

GMM=GAM( I-l, J,K) *GAM( I-l, J, K+1 ) 

1 /( ZCV( K ) *GAM( I-l , J , K+1 ) +ZCV( K+1 ) * 

2 GAM( I-l,J,K)+l.E-30) *XCVIP( I-l ) 

DIFF=YCV( J ) *2 . * ( GM+GMM) 

FLOW=YCV( J ) * ( FL+FLM) 

CALL PROFIL 

AKM( I , J , K+1 )=ACOF+AMAXl ( ZERO, FLOW) 

AKP ( I , J , K )=AKM( I , J , K+1 )-FLOW 
103 CONTINUE 

DO 104 J=2,M2 
DO 104 K=2,N2 

COEFFICIENTS AEAST AND AWEST 
AREA=YCV( J ) *ZCV( K ) 

FLOW=AREA*RHO( 1,J,K)*U(2,J,K) 

DI FF=AREA*GAM( 1 , J , K )/XCV( 2 ) 

CALL PROFIL 

AIM( 3 , J , K ) =ACOF+AMAXl ( ZERO, FLOW) 

FLOW=AREA*RHO( Ll,J,K)*U(Ll,J,K) 

DI FF=AREA*GAM( Ll , J , K )/XCV( L2 ) 

CALL PROFIL 

AIP( L2 , J , K )=ACOF+AMAXl ( ZERO, FLOW) -FLOW 

104 CONTINUE 

DO 105 1=3, L2 
DO 105 K=2,N2 

COEFFICIENTS ANORTH AND ASOUTH 

FL=XCVI ( I ) *V( I , 2 , K ) *RHO( I , 1 , K ) 

FLM=XCVIP( I-1)*V(I-1,2,K) *RHO( I-l , 1 , K) 

FLOW=ZCV( K) *( FL+FLM) 

GM=XCVI ( I ) *GAM( I , 1 , K ) +XCVIP( I-l ) *GAM{ I-l , 1 , K) 

DIFF=ZCV( K) *GM/YDIF( 2 ) 

CALL PROFIL 

AJM( I , 2 , K)=ACOF+AMAXl ( ZERO, FLOW ) 

FL=XCVI (I)*V(I,Ml,K) *RHO( I ,Ml , K ) 

FLM=XCVIP ( I-l ) *V( I-l , Ml , K ) *RHO( I-l , Ml , K ) 

FLOW=ZCV( K ) * ( FL+FLM) 

GM=XCVI ( I ) *GAM( I , Ml , K) +XCVIP( I-l ) *GAM( I-l ,Ml , K ) 

DIFF=ZCV( K) *GM/YDIF(Ml ) 

CALL PROFIL 

AJP( I , M2 , K )=ACOF+AMAXl (ZERO, FLOW) -FLOW 

105 CONTINUE 
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DO 106 1=3, L2 
DO 106 J=2,M2 
COEFFICIENTS AIN AND AOUT 

FL=XCVI (I)*W(I,J,2) *RHO( I , J , 1 ) 

FLM=XCVIP( I-1)*W(I-1,J,2) *RHO( I-l , J , 1 ) 

FLOW=YCV( J) *( FL+FLM) 

GM=XCVI ( I ) *GAN( I , J , 1 ) +XCVIP( I-l ) *GAM( I-l , J, 1 ) 

DIFF=YCV( J ) *GM/ZDIF( 2 ) 

CALL PROFIL 

AKM( I , J , 2 ) =ACOF+AMAXl (ZERO, FLOW) 

FL=XCVI ( I ) *W( I , J , N1 ) *RHO( I , J ,Nl ) 

FLM=XCVIP( I-l)*W(I-l,J,Nl) *RHO( I-l , J , Nl) 

FLOW=YCV( J ) * ( FL+FLM ) 

GM=XCVI ( I ) *GAN( I , J ,Nl )+XCVIP( I-l ) *GAM( I-l , J , Nl ) 

DIFF=YCV( J ) *GM/ZDIF(N1 ) 

CALL PROFIL 

AKP( I , J , N2 )=ACOF+AMAXl (ZERO, FLOW) -FLOW 
106 CONTINUE 

DO 107 1=3, L2 
DO 107 J=2,H2 
DO 107 K=2,N2 
VOL=YCV( J ) *XCVS( I ) *ZCV( K) 

APT=( RHO( I , J , K ) *XCVI ( I ) +RHO( I-l , J , K) *XCVIP( I-l ) ) 
l/( XCVS( I ) *DT) 

AP( I, J,K)=AP( I, J,K)-APT 

CON( I , J,K)=CON( I , J,K)+APT*U( I , J,K) 

AP( I , J,K)= 

1 ( -AP( I , J , K ) *VOL+AIP( I , J , K ) +AIH( I , J, K) +AJP( I , J , K ) +AJM( I , J , K) 

2 +AKP( I , J , K) +AKM( I , J , K) ) 

3/RELAX( NF ) 

CON( I , J,K)=CON( I , J,K) *VOL+REL*AP( I,J,K)*U(I,J,K) 

DU( I , J , K ) =VOL/XDIF( I ) 

DU( I , J , K) =DU( I , J , K)/AP( I , J , K ) 

107 CONTINUE 

DO 7099 1=1, LI 
DO 7099 J=1,M1 
DO 7099 K=l,Nl 
COFl( I,J,K,1)=AIP(I,J,K) 

COFl ( I , J , K , 2 ) =AIM( I , J , K ) 

COFl ( I , J , K , 3 ) =AJP( I , J , K ) 

COFl ( I , J , K , 4 ) =AJM( I , J , K ) 

COFl(I,J,K,5)=AKP(I,J,K) 

COFl ( I , J , K , 6 )=AKM( I , J , K ) 

COFl ( I , J , K, 7 )=AP( I , J , K) 

CONK I , J , K)=CON( I , J ,K) 

7099 CONTINUE 

C TEMPORARY USE OF PC(I,J) TO STORE UHAT 

DO 151 K=2,N2 
DO 151 J=2,M2 
DO 151 1=3, L2 

151 PC(I,J,K)=(AIP(I,J,K)*U(I+1,J,K)+AIM(I,J,K)*U(I-1,J,K) 

1 +AJP( I,J,K)*U(I,J+1,K) +AJM( I,J,K)*U(I,J-1,K) 

2 +AKP( I,J,K)*U(I,J,K+1) +AKM( I,J,K)*U(I,J,K-1) 

3 +CON( I , J , K) )/AP( I , J , K) 

100 CONTINUE 

C 

COEFFICIENTS FOR THE V EQUATION 

C 

CALL RESET 
NF=2 
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IF( .NOT.LSOLVE(NF) ) GO TO 200 

IST=2 

JST=3 

KST=2 

CALL DIFFUS 
REL=1 . -RELAX( NF ) 

DO 203 J=3,M2 
DO 203 1=2, L2 
DO 203 K=2,N2 

COEFFICIENTS ANORTH AND ASOUTH 

FL=V( I,J,K)*(FY(J) *RHO( I , J , K) +FYM( J ) *RHO( I , J-1 , K ) ) 

FLP=V( I,J+1,K)*(FY(J+1) *RHO( I,J+1,K)+FYM(J+1) *RHO( I , J , K ) ) 
FLOW=XCV( I)*ZCV(K)*0.5*( FL+FLP ) 

DIFF=XCV( I ) *ZCV( K ) *GAM( I , J , K )/YCV( J ) 

CALL PROFIL 

AJM( I , J+1 , K )=ACOF+AHAXl ( ZERO, FLOW) 

AJP( I , J , K )=AJM( I , J+1 , K )-FLOW 
COEFFICIENTS AEAST AND AWEST 

FL=YCVJ (J)*U(I+1,J,K)*(FX(I+1) *RHO( I+l , J , K ) +FXM( I+l)*RHO{I,J,K)) 
FLM=YCVJP ( J-1)*U(I+1,J-1,K)*(FX(I+1) *RHO( I+l , J-1 , K ) + FXM ( I+l ) * 

1 RHO( I , J-1 , K ) ) 

GM=GAM( I , J, K) *GAM( I+l , J, K) 

1 /( XCV( I ) *GAM( I+l , J , K )+XCV( I+l ) *GAM( I , J, K ) + 

2 1 . OE-30 ) *YCVJ( J ) 

GMM=GAM( I , J-1 ,K ) *GAM( I+l , J-1 , K ) 

1 /( XCV( I ) *GAM( I+l , J-1 , K ) +XCV( I+l ) * 

2 GAM( I,J-l,K)+l.E-30 ) *YCVJP( J-1 ) 

DIFF=ZCV( K ) *2 . * ( GM+GMM ) 

FLOW=ZCV( K ) * ( FL+FLM) 

CALL PROFIL 

AIM( I+l , J , K )=ACOF+AMAXl (ZERO, FLOW) 

AIP(I,J,K)=AIM( I+l, J,K) -FLOW 
COEFFICIENTS AIN AND AOUT 

FL=YCVJ( J )*W(I,J,K+1)*(FZ(K+1) *RHO( I,J,K+1)+FZM(K+1) *RHO( I , J, K ) ) 
FLM=YCVJP( J-1 ) *W( I , J-1 , K+1 ) * ( FZ( K+1 ) *RHO( I,J-1,K+1)+FZM(K+1)* 

1 RHO( I , J-1 , K ) ) 

GM=GAM( I , J, K ) *GAM( I , J, K+1 ) 

1 /( ZCV( K ) *GAM( I , J , K+1 ) +ZCV( K+1 ) *GAM( I , J , K) + 

2 1 . OE-30 ) *YCVJ( J ) 

GMM=GAM( I , J-1 , K ) *GAM( I , J-1 , K+1 ) 

1 /( ZCV( K ) *GAM( I , J-1 , K+1 ) +ZCV( K+1 ) * 

2 GAM( I,J-l,K)+l.E-30) *YCVJP( J-1 ) 

DIFF=XCV( I ) *2 . * (GM+GMM) 

FLOW=XCV( I ) * ( FL+FLM) 

CALL PROFIL 

AKM( I , J , K+1 )=ACOF+AMAXl ( ZERO, FLOW) 

AKP( I , J , K )=AKM( I , J , K+1 )-FLOW 
203 CONTINUE 

DO 204 1=2, L2 
DO 204 K=2,N2 

COEFFICIENTS ANORTH AND ASOUTH 
AREA=XCV( I ) *ZCV( K ) 

FLOW=AREA*RHO( I , 1 , K ) * V( I , 2 , K ) 

DIFF=AREA*GAM( 1,1, K )/YCV( 2 ) 

CALL PROFIL 

AJM( I, 3,K)=ACOF+AMAXl( ZERO, FLOW) 

FLOW=AREA*RHO( I , Ml , K ) * V( I , Ml , K ) 

DIFF=AREA*GAM( I , Ml , K ) /YCV ( M2 ) 

CALL PROFIL 

AJP( I ,M2 , K )=ACOF+AMAXl ( ZERO, FLOW) -FLOW 
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204 CONTINUE 

DO 205 J=3,n2 
DO 205 K=2,N2 

COEFFICIENTS AEAST AND AWEST 

FL=YCVJ ( J ) *U( 2 , J , K ) *RHO( 1 , J , K ) 

FLM=YCVJP( J-1 ) *U( 2 , J-1 , K ) *RHO( 1 , J-1 , K ) 

FLOW=ZCV( K) * ( FL+FLM) 

GM=YCVJ( J)*GAM(1, J,K)+YCVJP( J-1)*GAM(1, J-1,K) 

DIFF=ZCV( K) *GM/XDIF( 2 ) 

CALL PROFIL 

AIM( 2 , J , K ) =ACOF+AMAXl (ZERO, FLOW) 

FL=YCVJ (J)*U(Ll,J,K) *RHO( Ll , J , K ) 

FLM=YCVJP( J-1 )*U(L1,J-1,K) *RHO( Ll , J-1 , K) 

FLOW=ZCV( K ) * ( FL+FLM) 

GM=YCVJ ( J ) *GAM( Ll , J , K ) +YCVJP( J-1 ) *GAM( Ll , J-1 , K ) 

DIFF=ZCV( K) *GM/XDIF( Ll ) 

CALL PROFIL 

AI P( L2 , J , K )=ACOF+AMAXl ( ZERO, FLOW) -FLOW 

205 CONTINUE 

DO 206 1=2, L2 
DO 206 J=3,M2 
COEFFICIENTS AIN AND AOUT 

FL=YCVJ ( J ) *W( I , J , 2 ) *RHO{ I , J , 1 ) 

FLM=YCVJP( J-1 ) *W( I , J-1 , 2 ) *RHO( I , J-1 , 1 ) 

FLOW=XCV( I ) * ( FL+FLM) 

GM=YCVJ ( J ) *GAM( I , J , 1 ) +YCVJP{ J-1 ) *GAM( I , J-1 , 1 ) 

DIFF=XCV( I ) *GM/ZDI F( 2 ) 

CALL PROFIL 

AKM( I , J , 2 )=ACOF+AMAXl ( ZERO, FLOW) 

FL=YCVJ (J)*W(I,J,Nl) *RHO( I , J , N1 ) 

FLM=YCVJP( J-1 ) *W( I , J-1 ,N1 ) *RHO( I , J-1 , Nl ) 

FLOW=XCV( I ) * ( FL+FLM ) 

GM=YCVJ ( J ) *GAM( I , J ,Nl ) +YCVJP( J-1 ) *GAM( I , J-1 ,N1 ) 

DIFF=XCV( I ) *GM/ZDIF(Nl ) 

CALL PROFIL 

AKP( I , J , N2 ) =ACOF+AMAXl (ZERO, FLOW) -FLOW 

206 CONTINUE 

DO 207 1=2, L2 
DO 207 J=3,M2 
DO 207 K=2,N2 
VOL=XCV( I ) *YCVS( J ) *ZCV( K ) 

APT=( RHO( I , J , K ) * YCVJ ( J ) +RHO( I , J-1 , K) *YCVJP( J-1 ) ) 
l/( YCVS( J ) *DT) 

AP ( I , J , K ) =AP ( I , J , K ) -APT 

CON( I , J , K)=CON( I , J , K)+APT*V( I , J , K) 

AP( I , J,K)= 

1 ( -AP( I , J , K ) *VOL+AIP( I , J , K) +AIM( I , J , K) +AJP( I , J , K ) +AJH( I , J , K ) 

2 +AKP( I , J , K ) +AKM( I , J , K) ) 

3/RELAX( NF ) 

CON( I , J , K ) =CON( I , J , K ) *VOL+REL*AP( I,J,K)*V(I,J,K) 

DV( I , J , K )=VOL/YDIF( J ) 

DV( I , J , K ) =DV( I , J , K )/AP( I , J , K ) 

207 CONTINUE 

DO 8099 1=1, Ll 
DO 8099 J=1,M1 
DO 8099 K=1,N1 
COF2( I , J , K, 1 )=AIP( I , J, K) 

COF2 (I,J,K,2)=AIM(I,J,K) 

COF2( I , J , K, 3 )=AJP( I , J , K) 

COF2 ( I , J , K, 4 )=AJM( I , J , K) 
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C0F2 ( I , J , K , 5 ) =AKP ( I , J , K ) 

C0F2( I , j, K, 6 )=AKM( I, J,K) 

C0F2( I , J , K, 7 )=AP( I , J , K) 

C0N2( I , j , K)=C0N( I , j ,K) 

8099 CONTINUE 
200 CONTINUE 
C 

COEFFICIENTS FOR THE W EQUATION 

C 

CALL RESET 
NF=3 

IF( .NOT.LSOLVE(NF) ) GO TO 300 

IST=2 

JST=2 

KST=3 

CALL DIFFUS 
REL=1 .-RELAX(NF) 

DO 303 K=3,N2 
DO 303 J=2,M2 
DO 303 1=2, L2 

COEFFICIENTS AIN AND AOUT 

FL=W( I,J,K)*(FZ(K) *RHO( I , J , K ) +FZM ( K ) *RHO( I , J , K-1 ) ) 

FLP=W( I,J,K+1)*(FZ(K + 1) *RHO( I , J , K+1 ) +FZM( K + 1 ) *RHO( I , J , K ) ) 
FLOW=YCV( J ) *XCV( I ) *0 . 5* ( FL+FLP ) 

DIFF=YCV( J ) *XCV( I ) *GAM( I , J , K )/ZCV( K ) 

CALL PROFIL 

AKM( I , J , K+1 )=ACOF+AMAXl (ZERO, FLOW) 

AKP( I , J , K )=AKM( I , J , K+1 ) -FLOW 
COEFFICIENTS ANORTH AND ASOUTH 

FL=ZCVK( K )*V(I,J+1,K)*(FY(J+1) *RHO( I , J+1 , K ) +FYM( J+1 ) *RHO( I , J , K ) ) 
FLM=ZCVKP (K-1)*V(I,J+1,K-1)*(FY(J+1) *RHO( I , J+1 , K-1 ) +FYM( J+1 ) * 

1 RHO( I , J , K-1 ) ) 

GM=GAM( I, J,K)*GAM( I , J+1,K) 

1 /( YCV( J ) *GAM( I , J+1 , K)+YCV( J+1 ) *GAM( I , J , K) + 

2 1 . OE-30 ) *ZCVK( K) 

GMM=GAM( I , J , K-1 ) *GAM( I , J+1 , K-1 ) 

1 /( YCV( J ) *GAM( I , J+1 , K-1 ) +YCV( J+1 ) * 

2 GAM( I,J,K-l)+l.E-30) *ZCVKP( K-1 ) 

DIFF=XCV( I ) *2 . * ( GM+GMM) 

FLOW=XCV( I ) * ( FL+FLM) 

CALL PROFIL 

AJM( I , J+1 , K ) =ACOF+AMAXl ( ZERO, FLOW) 

AJP( I , J, K )=AJM( I , J+1 , K)-FLOW 
COEFFICIENTS AEAST AND AWEST 

FL=ZCVK( K )*U(I+1,J,K)*(FX(I+1) *RHO( I+l , J , K )+FXM( I+l ) *RHO( I , J , K) ) 
FLM=ZCVKP( K-1 )*U(I+1,J,K-1)*(FX(I+1) *RHO( I+l , J , K-1 )+FXM( I+l ) * 

1 RHO( I , J , K-1 ) ) 

GM=GAM( I , J , K) *GAM( I + l , J, K) 

1 /( XCV( I ) *GAM( I+l , J , K )+XCV( I+l ) *GAM( I , J , K )+ 

2 1 . OE-30 ) *ZCVK( K ) 

GMM=GAM( I , J , K-1 ) *GAM( I+l , J , K-1 ) 

1 /( XCV( I ) *GAM( I+l , J , K-1 )+XCV( I+l ) * 

2 GAM( I,J,K-l)+l.E-30) *ZCVKP( K-1 ) 

DIFF=YCV( J ) *2 . * (GM+GMM) 

FLOW=YCV( J ) * ( FL+FLM) 

CALL PROFIL 

AIM( I+l , J , K )=ACOF+AMAXl ( ZERO, FLOW) 

AIP(I,J,K)=AIM(I+1,J,K ) -FLOW 
303 CONTINUE 

DO 304 J=2,M2 
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DO 304 1=2, L2 

COEFFICIENTS AIN AND AOUT 
AREA=YCV( J ) *XCV( I ) 

FLOW=AREA*RHO( I,J,1)*W(I,J,2) 

DI FF=AREA*GAM( I , J , 1 )/ZCV( 2 ) 

CALL PROFIL 

AKM( I , J , 3 ) =ACOF+AMAXl ( ZERO, FLOW) 

FLOW=AREA*RHO( I , J , N1 ) *W( I , J , Nl) 

DIFF=AREA*GAM( I , J , Nl ) /ZCV ( N2 ) 

CALL PROFIL 

AKP( I , J ,N2 ) =ACOF+AMAXl ( ZERO, FLOW) -FLOW 

304 CONTINUE 

DO 305 1=2, L2 
DO 305 K=3,N2 

COEFFICIENTS ANORTH AND ASOUTH 

FL=ZCVK( K ) *V( I , 2 , K ) *RHO( I , 1 , K ) 

FLH=ZCVKP( K-1 ) *V( 1,2, K-1 ) *RHO( 1,1, K-1 ) 

FLOW=XCV( I ) * ( FL+FLM) 

GM=ZCVK( K ) *GAM( I , 1 , K ) +ZCVKP ( K-1 ) *GAM( I , 1 , K-1 ) 

DIFF=XCV( I ) *GM/YDIF( 2 ) 

CALL PROFIL 

AJM( I , 2 , K )=ACOF+AMAXl ( ZERO, FLOW) 

FL=ZCVK( K)*V(I,M1,K) *RHO( I , Ml , K ) 

FLM=ZCVKP( K-l)*V(I,Ml,K-l) *RHO( I , Ml , K-1 ) 

FLOW=XCV( I ) * ( FL+FLM) 

GM=ZCVK( K) *GAM( I ,Ml , K)+ZCVKP( K-1 ) *GAM( I ,Ml , K-1 ) 

DIFF = XCV( I ) *GM/YDI F( Ml ) 

CALL PROFIL 

AJP( I , M2, K)=ACOF+AMAXl( ZERO, FLOW) -FLOW 

305 CONTINUE 

DO 306 K=3,N2 
DO 306 J=2,M2 

COEFFICIENTS AEAST AND AWEST 

FL=ZCVK( K ) *U( 2 , J , K) *RHO( 1 , J , K) 

FLM = ZCVKP( K-1 ) *U ( 2 , J , K-1 ) *RHO( 1 , J , K-1 ) 

FLOW=YCV( J ) * ( FL+FLM) 

GM=ZCVK( K ) *GAM( 1 , J , K ) + ZCVKP ( K- 1 ) * GAM ( 1 , J , K-1 ) 

DIFF=YCV( J ) *GM/XDIF( 2 ) 

CALL PROFIL 

AIM( 2 , J , K ) =ACOF+AMAXl ( ZERO , FLOW ) 

FL=ZCVK( K)*U(L1,J,K) *RHO( Ll , J , K ) 

FLM=ZCVKP( K-1)*U(L1,J,K-1) *RHO( Ll , J , K-1 ) 

FLOW=YCV( J ) * ( FL+FLM) 

GM=ZCVK( K) *GAM( Ll , J , K )+ZCVKP( K-1 ) *GAM( Ll , J , K-1 ) 

DIFF=YCV( J ) *GM/XDIF( Ll ) 

CALL PROFIL 

AI P( L2 , J , K ) =ACOF+AMAXl ( ZERO, FLOW ) -FLOW 

306 CONTINUE 

DO 307 1=2, L2 
DO 307 J=2,M2 
DO 307 K=3,N2 
VOL=YCV( J ) *ZCVS( K) *XCV( I ) 

APT=( RHO( I , J , K) *ZCVK( K)+RHO( I , J , K-1 ) *ZCVKP( K-1 ) ) 

1/(ZCVS(K)*DT) 

AP( I , J ,K)=AP( I , J,K)-APT 

CON( I , J, K )=CON( I , J , K)+APT*W( I , J , K) 

AP( I , J , K ) = 

1 ( -AP ( I , J , K ) *VOL + AIP( I , J , K) +AIM( I , J , K ) fAJP( I , J , K) +AJM( I , J , K ) 

2 +AKP ( I , J , K ) +AKM( I , J , K ) ) 

3/RELAX( NF ) 
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CONI I , J, K)=CON( I , J, K) *VOL+REL*AP( I , J , K ) *W ( I , J , K ) 

DW( I, J,K)=VOL/ZDIF(K) 

DW( I , J , K) =DW( I , J , K)/AP( I , J , K ) 

307 CONTINUE 

DO 9099 1=1, LI 
DO 9099 J=1,M1 
DO 9099 K=1,N1 
COF3 { I , J , K, 1 )=AIP( I , J , K) 

COF3I I , J, K, 2 )=AIM( I , J,K) 

COF3I I, J,K, 3)=AJP( I,J,K) 

COF3 I I , J , K, 4 )=AJM{ I , J , K) 

COF3 ( I , J , K , 5 )=AKP( I , J , K ) 

COF3 I I , J , K, 6 )=AKM( I , J , K ) 

COF3(I,J,K,7) =AP( I , J , K) 

CON3 { I , J , K ) =CON ( I , J , K ) 

9099 CONTINUE 
300 CONTINUE 
C 

COEFFICIENTS FOR THE PRESSURE EQUATION 
C 

NF=NP 

IF( .NOT.LSOLVE(NF) ) GO TO 500 

IST=2 

JST=2 

KST=2 

DO 501 J=2,M2 
DO 501 K=2,N2 
AIM(2, J,K)=0.0 
AIPIL2, J,K)=0.0 

CONI 2 , J , K)=RHO| 1,J,K)*U|2,J,K) *YCV| J ) *ZCV| K ) 

501 CONTINUE 

DO 502 1=2, L2 
DO 502 K=2,N2 
AJMI I , 2 , K ) =0 . 0 
AJP| I , M2 , K ) =0 . 0 

CONI I , 2 , K ) =RHO| I,1,K)*V|I,2,K) *XCV{ I ) *ZCV| K ) 

502 CONTINUE 

DO 503 1=2, L2 
DO 503 J=2,M2 
AKMI I , J , 2 ) =0 . 0 
AKPI I , J , N2 ) =0 . 0 

CONI I , J , 2 )=RHO| I,J,1)*W|I,J,2) *XCV| I ) *YCV| J ) 

503 CONTINUE 
C 

DO 504 K=2,N2 
DO 504 J=2,m2 
DO 504 1=2, L2 
AREA=YCV| J ) *ZCV| K ) 

ARHO=AREA* | FX| 1 + 1 ) *RHO| I + 1 , J , K ) + FXM | I + 1 ) *RHO| I , J, K) ) 
FLOW=ARHO*PC| I+l, J,K) 

IF| I . EQ. L2 ) FLOW=ARHO*U| LI , J , K ) 

AIPI I , J, K) =ARHO*DU| I+l , J , K) 

AIM I I + l , J , K )=AIP| I, J,K) 

CONI I , J,K)=CON| I , J,K)-FLOW 
CONI I+l,J,K)=CON| I+l, J,K)+FLOW 
C 

AREA=XCV| I ) *ZCV| K) 

ARHO=AREA* | FY| J + 1 ) *RHO| I , J+1 , K ) +FYM| J + 1 ) * RHO I I , J , K ) ) 
VHAT = I COF2 I I , J + 1 , K, 1 ) *V| I + l , J + 1 , K ) 

1 +COF2I I,J+1,K,2)*V|I-1,J+1,K) 
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2 +C0F2 { I , J+1 , K , 3 ) *V( I , J + 2 , K ) 

3 +C0F2( I , J+1 , K, 4 ) *V( I , J, K) 

4 +C0F2 ( I , J+1 , K , 5 ) *V( I , J+1 , K+1 ) 

5 +C0F2 ( I , J+1 , K, 6 ) *V{ I , J+1 , K-1 ) 

6 +C0N2 ( I , J+1 , K ) )/( C0F2 (I,J+1,K,7)+1.0e-30) 
FLOW=ARHO*VHAT 

I F( J . EQ. M2 ) FLOW=ARHO*V( I , Ml , K) 

AJP{ I , J , K )=ARHO*DV( I , J+1 , K ) 

AJM{ I , J+1 ,K)=AJP{ I , J,K) 

CON( I , J,K)=CON( I ,J,K)-FLOW 
CON( I , J+1 , K)=CON( I , J+1 , K)+FLOW 

AREA=XCV( I ) *YCV( J ) 

ARHO=AREA* { FZ( K+1 ) *RHO( I , J , K+1 ) +FZM( K+1 ) *RHO( I , J , K) ) 
WHAT = ( C0F3 ( I , J, K+1 , 1 ) *W( I+l , J , K+1 ) 

1 +C0F3 ( I , J , K+1 , 2 ) *W( I-l , J , K+1 ) 

2 +C0F3 ( I , J , K+1 , 3 ) *W( I , J+1 , K+1 ) 

3 +C0F3 ( I , J , K+1 , 4 ) *W{ I , J-1 , K+1 ) 

4 +C0F3( I , J , K+1 , 5) *W( I , J , K+2 ) 

5 +C0F3 ( I , J , K+1 , 6 ) *W( I , J , K) 

6 +C0N3 ( I , J , K+1 ) )/(C0F3 ( I,J,K+l,7)+l.E-30) 
FLOW=ARHO*WHAT 

IF( K.EQ.N2)FL0W=ARH0*W( I , J,N1 ) 

AKP( I , J , K)=ARHO*DW{ I , J, K+1 ) 

AKM( I , J , K+1 )=AKP( I , J , K) 

CON( I , J, K)=CON( I , J,K)-FLOW 
CON{ I , J , K+1 )=FLOW 

AP( I , J , K)=AIP( I , J , K)+AIM( I , J , K ) +AJP( I , J , K ) +AJM( I , J , K) 

1 +AKP( I , J , K) +AKM( I , J , K ) 

504 CONTINUE 

DO 703 1=1 , Ll 
DO 703 J=1,M1 
DO 703 K=l,Nl 
COF4( I , J,K, 1)=AIP( I , J,K) 

COF4 (I,J,K,2)=AIM(I,J,K) 

COF4( I , J, K, 3 )=AJP( I , J, K) 

COF4 ( I , J , K, 4 )=AJM( I , J , K ) 

COF4( I , J,K,5)=AKP( I, J,K) 

COF4 ( I , J , K , 6 ) =AKM( I , J , K ) 

COF4( I , J,K,7 )=AP{ I, J,K) 

703 CONTINUE 

IF( ITER.LE. 1 ) GO TO 409 
DO 408 K=2,N2 
DO 408 J=2,M2 
DO 408 1=2, L2 

AP( I , J, K)=AP( I , J , K ) /RELAX (NP) 

CON( I , J , K )=CON( I , J , K )+( 1 . -RELAX ( NP ) ) *AP( I,J,K)*P(I,J,K) 

408 CONTINUE 

409 CONTINUE 
CALL TDMA 
NF=1 
IST=3 
JST=2 

KST = 2 

DO 704 1=1, Ll 
DO 704 J=1,M1 
DO 704 K=l,Nl 
AIP( I , J , K )=COFl ( I , J , K, 1 ) 

AIM{ I , J, K)=COFl{ I , J, K, 2 ) 

AJP( I , J , K )=COFl ( I , J , K, 3 ) 
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AJH( I , J , K) =C0F1 ( I , J , K, 4 ) 

AKP( I , J,K)=COFl( I , J,K, 5) 

AKM( I , J,K)=C0F1 ( I , J, K, 6 ) 

AP( I , J,K)=C0F1( I , J,K,7 ) 

CON( I , J,K)=C0N1( I , J,K) 

704 CONTINUE 

DO 413 K=2,N2 
DO 413 J=2,M2 
DO 413 1=3, L2 

413 CON( I , J , K)=CON( I , J, K)+DU( I , J , K) *AP( I,J,K)*(P(I-1,J,K)-P(I,J,K)) 
CALL TDMA 
NF = 2 
I ST=2 
JST=3 
KST=2 

DO 714 1=1, Ll 
DO 714 J=l,Ml 
DO 714 K=1,N1 
AIP( I , J , K )=COF2 ( I , J , K, 1 ) 

AIM( I, J,K)=C0F2( I , J,K,2) 

AJP ( I , J , K ) =COF2 U / J , K, 3 ) 

AJH(I,J,K)=COF2(I,J,K,4) 

AKP( I,J,K)=COF2(I,J,K,5) 

AKM( I , J,K)=C0F2( I , J,K,6) 

AP(I,J,K)=COF2(I,J,K,7) 

CON( I, J,K)=C0N2( I, J,K) 

714 CONTINUE 

DO 513 K=2,N2 
DO 513 1=2, L2 
DO 513 J=3,m2 

513 CON( I , J, K)=CON( I , J , K) +DV( I , J, K) *AP( I,J,K)*(P(I,J-1,K)-P(I,J,K)) 
CALL TDMA 
NF=3 
IST=2 
JST=2 
KST=3 

DO 724 1=1, Ll 
DO 724 J=l,Ml 
DO 724 K=l,Nl 
AIP(I,J,K)=COF3(I,J,K,l) 

AIM(I,J,K)=COF3(I,J,K,2) 

AJP( I , J , K ) =C0F3 ( I , J , K , 3 ) 

AJM( I , J , K ) =COF3( I , J , K, 4 ) 

AKP( I , J , K ) =COF3 ( I , J , K , 5 ) 

AKM( I , J , K) =COF3 ( I , J, K, 6 ) 

AP( I , J , K) =COF3( I , J , K, 7 ) 

CON( I , J , K )=CON3( I , J , K) 

724 CONTINUE 

DO 523 J=2,M2 
DO 523 1=2, l2 
DO 523 K=3,N2 

523 CON( I , J , K )=CON( I,J,K)+DW(I,J,K)*AP(I,J,K)*(P(I,J,K-1)-P(I,J,K)) 
CALL TDMA 
C 

COEFFICIENTS FOR THE PRESSURE CORRECTION EQUATION 

C 

DO 706 1=1, Ll 
DO 706 J=1,M1 
DO 706 K=1,N1 
AIP( I , J , K )=COF4 ( I , J , K , 1 ) 
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AIM( I , J , K ) =C0F4 ( I , J , K, 2 ) 

AJP( I , J , K)=C0F4 ( I , J , K, 3 ) 

AJM( I , J , K )=C0F4 ( I , J , K, 4 ) 

AKP( I , J , K )=C0F4 ( I , J , K , 5 ) 

AKM( I , J , K )=C0F4 ( I , J , K, 6 ) 

AP( I , J, K)=C0F4 ( I , J , K, 7 ) 

706 CONTINUE 
NF=4 

IF( .NOT.LSOLVE(NF) ) GO TO 500 

IST=2 

JST=2 

KST=2 

CALL DIFFUS 
SHAX=0 . 

SSUM=0 . 

DO 410 K=2,N2 
DO 410 J=2,M2 
DO 410 1=2, L2 
VOL=YCV( J ) *XCV( I ) *ZCV( K) 

410 CON( I , J, K)=CON( I , J , K) *VOL 
DO 701 K=2,n2 
DO 701 J=2,M2 

CON( 2 , J , K)=RHO( 1,J,K)*U(2,J,K) *YCV( J ) *ZCV( K) 

701 CONTINUE 

DO 702 1=2, L2 
DO 702 K=2,N2 

CON( I , 2 , K) =RHO( I,1,K)*V(I,2,K) *XCV( I ) *ZCV( K ) 

702 CONTINUE 

DO 153 1=2, L2 
DO 153 J=2,M2 

CON( I , J , 2 ) =RHO( I,J,1)*W(I,J,2) *XCV( I ) * YCV( J ) 

153 CONTINUE 
C 

DO 351 K=2,N2 
DO 351 J=2,M2 
DO 351 1 = 2, L2 
AREA=YCV( J ) *ZCV( K ) 

ARHO=AREA* ( FX( I + l ) *RHO( I + l , J , K ) + FXM( I + l ) *RHO( I , J , K ) ) 
FLOW=ARHO*U( I+l , J , K) 

CON( I , J , K)=CON( I , J , K) -FLOW 
CON( I+l , J , K )=CON( I+l , J , K) +FLOW 
C 

AREA=XCV( I ) *ZCV( K) 

ARHO=AREA* ( FY( J + 1 ) *RHO( I , J + 1 , K ) +FYM( J + 1 ) *RHO( I , J , K ) ) 
FLOW=ARHO*V( I , J+1 , K) 

CON( I , J , K)=CON( I , J , K ) -FLOW 
CON( I , J+1 , K) =CON( I , J+1 , K) +FLOW 
C 

AREA=XCV( I ) *YCV( J ) 

ARHO=AREA* ( FZ ( K+1 ) *RHO( I,J,K + 1)+FZM(K + 1) *RHO( I , J , K ) ) 
FLOW=ARHO*W( I , J , K+1 ) 

CON ( I , J , K ) =CON ( I , J , K ) -FLOW 
CON( I , J,K+1 )=FLOW 
PC( I , J , K ) = 0 . 

351 CONTINUE 

DO 352 1=2, L2 
DO 352 J=2,M2 
DO 352 K=2,n2 

SMAX=AMAX1 ( SMAX, ABS( CON( I , J , K) ) ) 

SSUM=SSUM+CON( I, J,K) 
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352 CONTINUE 
CALL TDHA 
C 

COME HERE TO CORRECT THE VELOCITIES 

C 

DO 511 K=2,N2 
DO 511 J=2,M2 
DO 511 1=2, L2 

IF(I.NE.2) U(I,J,K)=U(I,J,K)+DU(I,J,K)*(PC(I-1,J,K)-PC(I,J,K)) 
IF(J.NE.2) V(I,J,K)=V(I,J,K) +DV( I,J,K)*(PC(I,J-1,K)-PC(I,J,K)) 
IF(K.NE.2) W(I,J,K)=W(I,J,K)+DW(I,J,K)*(PC(I,J,K-1)-PC(I,J,K)) 
511 CONTINUE 
500 CONTINUE 
C 

COEFFICIENTS FOR OTHER EQUATIONS 

C 

CALL RESET 
IST=2 
JST = 2 
KST=2 

DO 600 NF=5,NFMAX 
IF( .NOT.LSOLVE(NF) ) GO TO 600 
CALL DIFFUS 
REL=1 .-RELAX(NF) 

DO 601 1=2, L2 
DO 601 J=2,M2 
DO 601 K=2,N2 

COEFFICIENTS AWEST AND AEAST 
AREA= YCV( J ) *ZCV( K ) 

FLOW=AREA*U( I+1,J,K)*(FX(I+1) *RHO( I+l , J , K ) +FXM( I+l ) *RHO( I , J , K ) ) 
DIFF=AREA*2 . *GAM( I , J , K ) *GAM( I+l , J , K )/( XCV( I ) *GAM( I+l , J , K) + 

1 XCV( I+l ) *GAM( I , J,K)+1 .OE-30) 

CALL PROFIL 

AIM( I+l , J , K )=ACOF+AMAXl ( ZERO, FLOW) 

AIP(I,J,K)=AIM(I+1,J, K )-FLOW 
COEFFICIENTS ANORTH AND ASOUTH 
AREA= XCV(I)*ZCV(K) 

FLOW=AREA*V( I,J + 1,K)*(FY(J + 1) *RHO( I , J+1 , K ) +FYM( J + 1 ) *RHO( I , J , K) ) 
DIFF=AREA*2 . *GAM( I , J, K ) *GAM( I , J+1 , K )/( YCV( J ) *GAM( I , J+1 , K ) + 

1 YCV( J+1) *GAM( I , J , K) +1 . OE-30 ) 

CALL PROFIL 

AJM( I , J+1 , K) =ACOF+AMAXl ( ZERO, FLOW) 

AJP( I , J, K )=AJM( I , J+1 , K )-FLOW 
COEFFICIENTS AOUT AND AINTO 
AREA= YCV( J)*XCV( I ) 

FLOW=AREA*W( I , J , K+1 ) * ( FZ ( K+1 ) *RHO( I , J , K+1 )+FZM( K+1 ) *RHO( I , J , K) ) 
DIFF=AREA*2 , *GAM( I , J , K ) *GAM( I , J, K+1 )/( ZCV( K) *GAM( I , J , K+1 )+ 

1 ZCV( K+1 ) *GAM( I , J , K)+l . 0E-30 ) 

CALL PROFIL 

AKM( I , J , K+1 ) =ACOF+AMAXl ( ZERO, FLOW) 

AKP( I , J , K )=AKM( I , J , K+1 )-FLOW 
601 CONTINUE 

DO 610 J=2,M2 
DO 610 K=2,N2 

COEFFICIENTS AWEST AND AEAST 
AREA = YCV( J ) *ZCV( K ) 

FLOW=AREA*U( 2 , J , K ) *RHO( 1 , J, K) 

DI FF=AREA*GAM( 1 , J , K )/XDI F( 2 ) 

CALL PROFIL 

AIM( 2 , J , K ) =ACOF+AMAXl ( ZERO, FLOW) 
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FLOW=AREA*U( Ll , J , K) *RHO( Ll , J , K) 

DIFF=AREA*GAM( Ll , J , K)/XDIF( Ll ) 

CALL PROFIL 

AIP ( L2 , J , K ) =ACOF+AMAXl ( ZERO, FLOW) -FLOW 

610 CONTINUE 

DO 611 1=2, L2 
DO 611 K=2,N2 

COEFFICIENTS ANORTH AND ASOUTH 
AREA= XCV( I ) *ZCV( K) 

FLOW=AREA*V( I , 2, K) *RHO( I , 1 , K) 

DI FF=AREA*GAM( I,I,K)/YDIF(2) 

CALL PROFIL 

AJM( 1,2, K )=ACOF+AHAXl ( ZERO, FLOW) 

FLOW=AREA* V( I , Ml , K ) *RHO( I ,Ml , K ) 

DI FF=AREA*GAN( I,Ml,K)/YDIF(Ml) 

CALL PROFIL 

AJP( I , M2 , K )=ACOF+AMAXl ( ZERO, FLOW) -FLOW 

611 CONTINUE 

DO 612 1=2, L2 
DO 612 J=2,M2 

COEFFICIENTS AOUT AND AINTO 
AREA= YCV( J ) *XCV( I ) 

FLOW=AREA*W( I , J , 2 ) *RHO( I , J , 1 ) 

DIFF=AREA*GAM( I,J,1)/ZDIF(2) 

CALL PROFIL 

AKM( I , J, 2 )=ACOF+AMAXl (ZERO, FLOW) 

FLOW=AREA*W( I , J , Nl ) *RHO( I , J , Nl ) 

DIFF=AREA*GAM( I,J,Nl)/ZDIF(Nl) 

CALL PROFIL 

AKP( I , J ,N2 )=ACOF+AMAXl ( ZERO, FLOW) -FLOW 

612 CONTINUE 

DO 3987 1=1, Ll 
DO 3987 J=l,Ml 
DO 3987 K=l,Nl 
VOL=YCV( J ) *XCV( I ) *ZCV( K ) 

APT=RHO( I , J , K )/DT 

AP( I, J,K)=AP( I , J,K)-APT 

CON( I , J , K )=CON( I , J , K) +APT*F( I , J , K, NF) 

AP( I , J, K) 

1 =( -AP( I , J , K ) *VOL+AIP( I , J , K )+AIM( I , J , K ) +AJP( I , J , K ) +AJM( I , J , K ) 

2 +AKM( 1 , J , K )+AKP( I , J , K) ) 

3/RELAX( NF ) 

CON( I , J , K )=CON( I , J , K) *VOL+REL*AP( I,J,K)*F(I,J,K,NF) 

C PRINT*,I,J,K 

C PRINT* , AIM ( I , J , K ) , AJM( I , J , K ) , AKM( I , J , K) 

C PRINT*,AIP(I,J,K),AJP(I,J,K) ,AKP( I, J,K) 

C PRINT* , AP( I , J , K ) , CON( I , J , K ) 

3987 CONTINUE 
CALL TDMA 
600 CONTINUE 

TIME=TIME+DT 

ITER=ITER+1 

IF( ITER. GE. LAST) LSTOP=.TRUE. 

RETURN 

END 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
SUBROUTINE OPTION 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
LOGICAL LSOLVE, LPRINT, LBLK, LSTOP 

COMMON F(17,17,13,5),P(17,17,13), RHO( 17,17,13) ,GAM( 17 , 17 , 13 ) , 
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n n n n 



1 CON( 17,17,13),AKP(17,17,13) , AKM{ 1 7 , 1 7 , 1 3 ) ,AP(17,17,13), 

2 AIP(17,17,13),AIM(17,17,13),AJP(17,17,13) , AJM( 17,17,13) 

COMMON 

3 X(17),XU(17),XDIF(17)> XCV( 17 ) ,XCVS( 17 ) , 

4 Y(17),YV(17),YDIF(17) , YCV( 17) , YCVS( 17 ) , 

5 Z(17),ZW(17),ZDIF(17), ZCV( 17 ) , ZCVS( 17 ) , 

6 YCVR( 17 ) , YCVRS( 17 ) , ARX( 17 ) ,ARXJ( 17 ) , ARXJP( 17 ) , 

7 R( 17 ) , RMN( 17 ) , SX( 17 ) ,SXMN( 17 ) ,XCVI ( 17 ) ,XCVIP( 17 ) , 

0 YCVJ( 17 ) , YCVJP( 17 ) , ZCVK( 17 ) ,ZCVKP( 17 ) 

COMMON DU (17, 17, 13), DV (17, 17, 13), DW (17, 17, 13), FV (17), FVP( 17 ) , 

1 FX( 17 ) , FXM( 17 ) , FY( 17 ) , FYM( 1 7 ) , PT ( 1 7 ) , QT ( 1 7 ) , 

2 FZ U7 ) , FZM( 17 ) 

COMMON/INDX/NF,NFMAX, NP, NRHO, NGAM, L1,L2,L3,M1,M2,M3,N1,N2,N3, 

1 1ST, JST, KST, ITER, LAST, TITLE ( 13 ) , RELAX ( 13 ) , TI ME , DT , XL , YL , Z L , 

2l PREF , JPREF , KPREF , LSOLVE( 10), LPRINT( 1 3 ) , LBLK ( 1 1 ) , MODE, NTIMES ( 10 ) , 
3RHOCON, ZERO 

COMMON/INDX/RELAX( 13),LPRINT(13),LBLK(11), NTIMES ( 10 ) , 

1LS0LVE( 10 ) , TIME, DT, XL, YL, ZL, S , RHOCON, ZERO, 
2NF,NFMAX,NP,NRH0,NGAM,L1,L2,L3,M1,M2,M3,N1,N2,N3, 

3IST, JST,KST, ITER, LAST, 

4 I PREF , JPREF , KPREF , MODE 
COMMON/HEAD I N/TI TLE 
CHARACTER*10 TITLE(13) 

COMMON/CNTL/LSTOP 
COMMON/SORC/SMAX, SSUM 
COMMON/COEF/FLOW, DI FF , ACOF 

DIMENSION U(17,17,13),V(17,17,13),W(17,17,13),PC(17,17,13) 
DIMENSION T(17,17,13) 

EQUIVALENCE( F( 1,1,1,1),U(1,1,1)),(F(1,1,1,2),V(1,1,1)), 

1 (F( 1, 1, 1, 3 ) ,W( 1,1,1) ),(F( 1,1,1, 4 ),PC( 1,1,1 ) ) 

2 , (F(l, 1,1,5), T(l, 1,1) ) 



10 FORMAT ( 26(1H*),3X,A10,3X,26(1H*)) 

20 FORMAT(lX,4H I =,I6,6I9) 

3 0 FORMAT( lx , IHJ ) 

40 FORMAT( IX, 12 , 3X, IP, 7E9 . 2 ) 

50 FORMAT UH ) 

51 FORMAT(lX,'I =' , 2X, 7 ( 14 , 5X) ) 

52 FORMAT ( IX, 'X =',1P,7E9.2) 

53 FORMAT('TH =',1P,7E9.2) 

54 FORMAT ( lx, 'J =' , 2X, 7 ( 14 , 5X) ) 

55 FORMAT(lX,'Y =',1P,7E9.2) 

56 FORMAT(lX,'K = ' , 2X , 7 ( I 4 , 5X ) ) 

57 FORMAT(lX,'Z =',1P,7E9.2) 

59 FORMAT(lX,'K =',2X,l4) 

ENTRY UHESH 

(2*Ai»c;^Ai*t**T»ti»t>:*****AAik*****7k**ik*T*t**7k7k*AAi»tT»t*7kAA**** 

XU( 2 )=0 . 

DX=XL/FLOAT( Ll-2 ) 

DO 1 1=3, LI 

1 XU( I ) =XU( I-l ) +DX 
YV( 2)=0. 

DY=YL/FLOAT( Ml-2 ) 

DO 2 J=3,H1 

2 YV( J )=YV( J-1 ) +DY 
ZW( 2 ) =0 . 

DZ=ZL/FLOAT( N1-2 ) 

DO 3 K=3,N1 

3 ZW( K)=ZW(K-1 )+DZ 



RETURN 

(^*AAA*AAA***A********A***AAAAAA*AAAAA*AAAA***AAAA 

ENTRY PRINT 

qaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa 
IF( .NOT. LPRINT( 3 ) ) GO TO 80 
80 CONTINUE 



IF( .NOT.LPRINT(NP) ) GO TO 90 
C 

CONSTRUCT BOUNDARY PRESSURES BY EXTRAPOLATION 
C 

DO 91 K=2,N2 
DO 91 J=2,m2 

P(1,J,K)=(P(2,J,K) *XCVS (3)-P(3,J,K)*XDIF(2) )/XDI F( 3 ) 

91 P(Ll,J,K)=(P(L2,J,K) *XCVS (L2)-P(L3,J,K)*XDIF(L1))/XDIF(L2) 
DO 92 K=2,N2 

DO 92 1=2, L2 

P(I,1,K)=(P(I,2,K) *YCVS( 3)-P(I,3,K)*YDIF(2) )/YDIF( 3 ) 

92 P(I,Ml,K)=(P(I,M2,K) * YCVS( M2 ) -P ( I , M3 , K ) * YDI F ( Ml ) )/YDIF( H2 ) 
DO 93 J=2,M2 

DO 93 1=2, L2 

P(I,J,1)=(P(I,J,2) *ZCVS( 3)-P(I,J,3)*ZDIF(2))/ZDIF(3) 

93 P(I,J,N1)=(P(I,J,N2) *ZCVS( N2 )-P( I,J,N3)*ZDIF(Nl) )/ZDIF(N2) 
DO 94 K=2,N2 

P(1,1,K)=P(2,1,K)+P(1,2,K)-P(2,2,K) 

94 CONTINUE 

DO 95 J=2,M2 

P(1,J,1)=P(2,J,1)+P(1,J,2)-P(2,J,2) 

95 CONTINUE 

DO 96 1=2, L2 

P( I , 1 , 1 )=P( I , 2, 1 )+P( I , 1 , 2)-P( I , 2, 2) 

96 CONTINUE 

P(l,l,l)=(P(l,l,2)+P(l,2,l)+P(2,l,l))/3.0 

P(Ll,l,l)=(P(L2,l,l)+P(Ll,2,l)+P(Ll,l,2))/3.0 

P(l,l,Nl)=(P(l,l,N2)+P(l,2,Nl)+P(2,l,Nl))/3.0 

P(l,Ml,l)=(P(l,M2,l)+P(l,Ml,2)+P(2,Ml,l))/3.0 

P(L1,H1,1)=(P(L2,M1,1)+P(L1,m2,1)+P(L1,m1,2))/3.0 

P(1,M1,N1)=(P(2,M1,N1)+P(1,M2,N1)+P(1,M1,N2))/3.0 

P( LI , Ml , Nl ) =( P( L2 ,M1 , N1 ) +P( LI , M2 , N1 ) +P( LI , Ml , N2 ) )/3 . 0 

PREF=P(IPREF,JPREF,KPREF) 

DO 97 K=l,Nl 
DO 97 J=l,Ml 
DO 97 1=1, LI 

97 P( I , J,K)=P( I , J,K)-PREF 
90 CONTINUE 

C 

PRINT 50 
IEND=0 

301 IF( lEND.EQ.Ll ) GO TO 310 
IBEG=IEND+1 
IEND=IEND+7 
IEND=MIN0( lEND, Ll ) 

PRINT 50 

PRINT 51 , ( I , I=IBEG, lEND) 

I F( MODE . EQ . 3 ) GO TO 302 
PRINT 52, (X( I ) , I=IBEG, lEND) 

GO TO 303 

302 PRINT 53 , ( X( I ) , I=IBEG, lEND) 

303 GO TO 301 
310 JEND=0 
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PRINT 50 

311 IF( JEND. EQ . Ml ) GO TO 320 
JBEG=JEND+1 
JEND=JEND+7 
JEND=MIN0 ( JEND, Ml ) 

PRINT 50 

PRINT 54 , ( J, J = JBEG, JEND) 

PRINT 55, (Y( J) , J=JBEG, JEND) 

GO TO 311 

320 KEND=0 
PRINT 50 

321 IF( KEND. EQ. Nl ) GO TO 330 
KBEG=KEND+1 
KEND=KEND+7 

KEND=MIN0 ( KEND, Nl ) 

PRINT 50 

PRINT 56 , ( K, K=KBEG, KEND) 

PRINT 57, (Z(K) ,K=KBEG,KEND) 

GO TO 321 
330 CONTINUE 
C 

DO 999 NF=1,NGAM 

IF( .NOT.LPRINT(NF) ) GO TO 999 

PRINT 50 

PRINT 10,TITLE(NF) 

DO 998 K=l,Nl 
PRINT 50 
PRINT 59, K 
IFST=1 
JFST=1 
KFST=1 

IF(NF. EQ. 1 .OR .NF. EQ. 4 ) IFST=2 
IF( NF. EQ. 2 .OR.NF. EQ. 4 ) JFST=2 
IF( NF. EQ. 3 .OR.NF. EQ. 4 ) KFST=2 
IBEG=IFST-7 
110 CONTINUE 

IBEG=IBEG+7 
I END=IBEG+6 
IEND=MIN0( IEND,Ll ) 

PRINT 50 

PRINT 20 , ( I , I = IBEG, lEND) 

PRINT 30 
JFL=JFST+Ml 
DO 115 JJ=JFST,Ml 
J=JFL-JJ 

PRINT 40,J,(F(I,J,K,NF),I = IBEG, I END) 

115 CONTINUE 

IF( lEND.LT.Ll ) GO TO 110 

998 CONTINUE 

999 CONTINUE 
RETURN 
END 

C 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

C PROBLEM DEPENDENT PORTION C 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

BLOCK DATA 

LOGICAL LSOLVE, LPRINT, LBLK, LSTOP 

COMMON F(17,17,13,5),P(17,17,13) ,RHO( 17, 17,13) , GAM ( 1 7 , 1 7 , 1 3 ) , 

1 CON (17, 17, 13) ,AKP(17,17,13) , AKM( 17 , 17 , 13 ) ,AP(17,17,13), 



40 



2 AIP(17,17,13),AIM(17,17,13),AJP(17,17,13) , AJM ( 1 7 , 1 7 , 1 3 ) 
COMMON 

3 X{17),XU(17),XDIF(17) ,XCV( 17 ) ,XCVS( 17 ) , 

4 Y(17),YV(17),YDIF(17) , YCV( 17 ) , YCVS( 17 ) , 

5 Z(17),ZW(17),ZDIF(17),ZCV(17), ZCVS( 17 ) , 

6 YCVR( 17 ) , YCVRS( 17 ) , ARX( 17 ) , ARX J ( 17 ) , ARXJP( 17 ) , 

7 R( 17 ) , RMN( 17 ) , SX( 17 ) ,SXMN( 17 ) ,XCVI ( 17 ) ,XCVIP( 17 ) , 

8 YCVJ(17) ,YCVJP(17) ,ZCVK(17) ,ZCVKP(17) 

COMMON DU( 17, 17 , 13 ) ,DV( 17, 17 , 13) ,DW( 17, 17, 13) ,FV( 17 ) , FVP( 17 ) , 

1 FX( 17 ) , FXM( 17 ) , FY( 17 ) , FYM( 1 7 ) , PT ( 1 7 ) , QT ( 1 7 ) , 

2 FZ( 17 ) , FZM( 17 ) 

COMMON/I NDX/RELAX( 13 ) , LPRINT( 13 ) , LBLK( 1 1 ) , NTIMES ( 1 0 ) , 

1LS0LVE( 10 ) , TIME , DT,XL, YL, ZL, S , RHOCON, ZERO, 

2NF,NFMAX, NP,NRHO,NGAM, LI , L2 , L3 , Ml , M2 , M3 , Nl , N2 , N3 , 
3IST,JST,KST, ITER, LAST, 

4IPREF,JPREF, KPREF , MODE 
COMMON/HEAD IN/TITLE 
CHARACTER*10 TITLE(13) 

COMMON/CNTL/LSTOP 
COMMON/SORC/SMAX, SSUM 
COMMON/COEF/FLOW , DI FF , ACOF 

DIMENSION U(17,17,13),V(17,17,13),W(17,17,13),PC(17,17,13) 
DIMENSION T( 17 , 17 , 13 ) 

EQUIVALENCE(F(1,1,1,1),U(1,1,1)),(F(1,1,1,2),V( 1,1,1)), 

1 (F(1,1,1,3),W(1,1,1)),(F(1,1,1,4),PC(1,1,D) 

2 ,(F(1,1,1,5),T(1,1,D) 

DATA NFMAX,NP,NRH0,NGAM/5,6,7,8/ 

DATA LSTOP, LSOLVE, LPRINT/2 4 * . FALSE./ 

DATA MODE , TIME, ITER/1 , 0 ., 0/ 

DATA RELAX, NTIMES/13*! . 0 , 10*1/ 

DATA LBLK/1 1 * . TRUE ./ 

DATA DT, IPREF, JPREF, KPREF, RHOCON/1 . E+ 1 0 , 1 , 1 , 1 , 1 . 0/ 



C NF=1,2,3 STAND FOR U,V AND W VELOCITIES. 

2 NF=6 IS FOR PRESSURE 
2 NF=5 IS FOR TEMPERATURE 

2 LSOLVE=TRUE SOLVES THAT PARTICULAR PHI 
DATA ( LSOLVE{ I ) , 1=1 , 6 )/6* . TRUE ./ 

2 LPRINT( NF ) =TRUE PRINTS VARIABLE ASSOCIATED WITH NF ON CALLING PRINT 
DATA LPRINT( 5 )/l * . TRUE ./ 

2 TERMINATE ITERATIONS AT ITER=LAST 
DATA LAST/100/ 

2 UNDERELAXATION FACTORS 

DATA RELAX( 1 ) , RELAX( 2 ) , RELAX( 3 ) , RELAX( 4 )/0 .5, 0.5, 0.5, 1.0/ 

DATA RELAX( 5 ) ,RELAX( 6 )/l . 0 , 1 . 0/ 

P TITLES FOR THE FIELD PRINTOUTS 

DATA TITLE( 1 ) , TITLE( 2 ) , TITLE( 3 ) , TITLE{ 5 ) , TITLE( 6 )/ 

I'U' , 'V' , 'W' , 'T' , 'P'/ 

2 NUMBER OF SWEEPS IN THE LINE-BY-LINE TDMA ALGORITHM 
DATA NTIMES( 3) ,NTIMES( 4 )/2*5/ 

DATA ZERO/0./ 

END 






SUBROUTINE USE 

2 IF LARGER NUMBER OF GRID POINTS IS TO BE USED, THE DIMENSION 
2 STATEMENTS MUST BE CHANGED THROUGH THE PROGRAM TO ACCOMODATE 
2 VALUES GRAETER THAN (17,17,13) ETC. 

LOGICAL LREAD,LWRITE 



I 
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no on 



LOGICAL LSOLVE, LPRINT, LBLK, LSTOP 

COMMON F(17,17,13,5),P(17,17,13), RM0( 17 , 17,13) ,GAM( 17 ,17,13), 

1 CON( 17,17,13) ,AKP(17,17,13),AKM(17,17,13) ,AP(17,17,13), 

2 AIP(17,17,13),AIM(17,17,13) ,AJP( 17,17,13),AJM(17,17,13) 

COMMON 

3 X(17),XU(17),XDIF(17) ,XCV( 17 ) , XCVS ( 17 ) , 

4 Y(17),YV(17),YDIF(17) , YCV( 17 ) , YCVS( 17 ) , 

5 Z( 17 ) ,ZW( 17 ) ,ZDIF( 17 ) ,ZCV( 17 ) ,ZCVS( 17 ) , 

6 YCVR( 17 ) , YCVRS( 17 ) ,ARX( 17 ) , ARXJ( 17 ) , ARXJP( 17 ) , 

7 R( 17 ) , RMN( 17 ) ,SX( 17 ) , SXMN( 17 ) ,XCVI ( 17 ) ,XCVIP( 17 ) , 

8 YCVJ( 17 ) , YCVJP( 17 ) , ZCVK( 17 ) , ZCVKP( 17 ) 

COMMON DU(17,17,13),DV(17,17,13),DW(17,17,13),FV(17), FVP( 17 ) , 

1 FX( 17 ) , FXM( 17 ) , FY( 17 ) , FYM( 17 ) , PT( 17 ) , QT( 17 ) , 

2 FZ( 17 ) , FZM( 17 ) 

COMMON/INDX/RELAX( 13 ) , LPRINT( 13 ) , LBLK( 11 ) , NTIMES ( 10 ) , 

1LS0LVE( 10) ,TIME,DT,XL,YL,ZL,S, RHOCON , ZERO, 

2NF, NFMAX,NP, NRHO,NGAM, LI , L2 , L3 , Ml , M2 , M3 , Nl , N2 , N3 , 

3IST,JST,KST, ITER, LAST, 

4IPREF, JPREF,KPREF,MODE 
COMMON/HEAD IN/TITLE 
CHARACTER*10 TITLE(13) 

COMMON/CNTL/LSTOP 
COMMON/SORC/SMAX , SSUM 
COMMON/COEF/FLOW,DIFF , ACOF 

DIMENSION U(17,17,13),V(17,17,13),W(17,17,13),PC(17,17,13) 
DIMENSION T(17,17,13) , BUOY( 17 ,17,13) 

C SPECIFY EQUIVALENCE FOR QUANTITIES INSIDE AND OUTSIDE SUBROUTINE USE 
C T=TEMPERATURE, PC=PRESSURE CORRECTION 

EQUIVALENCE( F( 1 ,1,1,1),U(1,1,1)),(F(1,1,1,2),V(1,1,1)), 

1 (F(1,1,1,3),W(1,1,1)),(F(1,1,1,4),PC(1,1,D) 

2 , (F(l, 1,1,5), T(l, 1,1) ) 

DIMENSION ANUZ ( 13 ) ,ANU3PT( 13 ) 

CHARACTER*! AREAD 

'k-k-k'k'k'kif'k'k-k'k-k'k-k'k-k'k'k-k-k'kif-k'kif-k'k'k-kif-k'k-A'kif'k'k-k-k'k'k-kif'k'k'kif 

ENTRY MESH 

DOMAIN LENGHTS IN THE 3 DIRECTIONS 
XL=1 . 0 
YL=1 . 0 
ZL=1 . 0 

C NUMBER OF GRID POINTS IN THE 3 DIRECTIONS 
Ll = 17 
Ml = 17 
Nl = 13 

C INVOKE THE UNIFORM MESH OPTION 
CALL UMESH 
RETURN 
C 

ENTRY BEGIN 

C ITERATIONS STOP AT ITER=LAST 

WRITE( * , * ) 'NO. ITERATIONS' 

READ( * , * ) LAST 
C READ RAYLEIGH NUMBER 
WRITE( * , * ) ' RALI= ' 

READ( * , * )RA 
C READ PRNDTL NUMBER 

WRITE( * , * ) ' PRANTL= ' 

READ( * , * ) PR 

C DETERMINE IF WANT TO USE A PREVIOUSLY COMPUTED SOLUTION AS 



42 



n o 



C AN INITIAL GUESS 

WRITE( * , * ) ' READ FROM INPUT FILE (0/1)' 

READ( * , * ) IREAD 
I F( IREAD . EQ. 1 ) LREAD= . TRUE . 

IF{ IREAD. EQ. 0 ) LREAD= . FALSE. 

C PROVIDE INITIAL GUESS. THE PROGRAM SOLVES FOR THE INTERIOR POINTS 
C ONLY. HENCE THE BOUNDARY VALUSE FOR THE TEMPERATURES AT THE HOT 
C AND COLD BOUNDARIES HAVE BEEN ALREADY SPECIFIED. 

DO 101 1=1, Ll 
DO 101 J=1,M1 
DO 101 K=l,Nl 
U( I , J, K)=0 . 

V( I , J, K)=0 . 

W( I , J , K ) =0 . 

T( I, J,K)= 0. 

C HOT WALL AT X=0. 

T(1,J,K)=1.0 
101 CONTINUE 

I F( . NOT . LREAD ) RETURN 
C READ DATA FROM INPUT FILE 
REWIND (4) 

READ ( 4 ) X , Y , Z , XU , YV , ZW , U , V , W , P , T 
RETURN 

ENTRY VARRHO 

(I^AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

RETURN 

C 

^aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa 
C INCORPORATE BOUNDARY CONDITIONS 
ENTRY BNDRY 
DO 864 1=1, Ll 
DO 864 J=l,Ml 

C UPDATE VELOCITIES AND TEMPERATURE AT THE SYMMETRY LINE (Z=ZL) 

C ACORDING TO THE ZERO GRADIENT BOUNDARY CONDITION FOR U,V AND T 
U( I , J , Nl )=U(I , J ,N2 ) 

V( I , J ,Nl )=V( I , J ,N2 ) 

T( I , J ,N1 )=T( I , J ,N2 ) 

C ADIABATIC SIDE (Z=0) 

T(I,J,1)=T(I,J,2) 

864 CONTINUE 

DO 865 1=1, Ll 
DO 865 K=1,N1 
C ADIABATIC BOTTOM ( Y=0 ) 

T(I,1,K)=T(I,2,K) 

C ADIABATIC TOP (Y=YL) 

T(I ,M1,K)=T( I ,M2,K) 

865 CONTINUE 

RETURN 

C 

ENTRY PRTOUT 
IF{ ITER.NE . 0 ) GO TO 400 
PRINT 401 

401 FORMAT { IX , ' S I MPLER ' , // ) 

400 CONTINUE 

C MONITOR SSUM,SMAX AND OTHER QUANTITES AS ITERATIONS PROCEED 
C ON CONVERGENCE, SMAX SHOULD BE VERY SMALL (LESS THAN l.OE-04) 

C SSUM SHOLD ACHIEVE A SMALL VALUE WITHIN A FEW ITERATIONS, WELL 

C BEFORE CONVERGENCE. SSUM WILL NOT BE SMALL IF THE BOUNDARY 
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C CONDITIONS ARE NOT WRITTEN CORRECTLY. 

PRINT 403 , ITER, SSUM, SMAX, U( 5,5,9),V(5,5,5),T(8,6,5) 

403 FORMAT (I6,5E12.4) 

COMPUTE LOCAL AND GLOBAL NUSSELT NUMBERS AS ITERATIONS PROCEED 
IF( ( ITER. EQ. 250 ) 

2 .OR. ( ITER. EQ. 500 ) 

3 .OR. ( ITER. EQ. 750 ) 

4 .OR. ( ITER. EQ. 1000 ) 

5 . OR . ( ITER . EQ . 1250 ) 

6 . OR. ( ITER . EQ. 1500 ) 

7 . OR . ( ITER . EQ . 1750 ) 

8 . OR . ( ITER . EQ . 2000 ) 

1 .OR. ( ITER. EQ. 2250 ) 

2 .OR. ( ITER. EQ. 2500 ) 

3 .OR . ( ITER . EQ . 2750 ) 

4 .OR. ( ITER. EQ. 3000 ) 

5 . OR , ( ITER . EQ , LAST ) ) THEN 

WRITE ( 9 , * ) ' ITER, SSUM, SMAX, U( 5,5,9),V(5,5,5),T(8,6,5)' 

WRITE (9,*)ITER,SSUM, SMAX,U( 5,5,9),V(5,5,5),T(8,6,5) 

COMPUTE AVERAGE NUSSELT NUMBER 
ANUHOT=0 . 0 
ANUCLD=0 . 0 
DO 666 K=l,Nl 
ANUZ( K)=0 . 0 
DO 667 J=1,M1 

COMPUTE AVERAGE NUSSELT NUMBER AT THE HOT AND COLD WALL 

ANUHOT=ANUHOT+( T( 1,J,K)-T(2,J,K) ) *YCV( J )*ZCV(K)/(XDIF(2)*YL*ZL) 
ANUCLD=ANUCLD+ (T(M2,J,K)-T(Ml,J,K) ) *YCV( J )*ZCV(K)/(XDIF(Ll)*YL*ZL) 
COMPUTE LOCAL NUSSELT NUMBER AT THE HOT WALL 

ANUZ( K )=ANUZ (K)+(T(1,J,K)-T(2,J,K)) *YCV( J )/( XDIF( 2 ) *YL) 

1 *YCV( J)/(XDIF( 2) *YL) 

667 CONTINUE 

C WRITE NUSSELT NUMBERS TO AN OUTPUT FILE 

WRITE(9,*)'Z(K)=',Z(K), ' ANUZ( K) =' , ANUZ( K) 

666 CONTINUE 

WRITE( 9 , * ) ' ANUHOT =',ANUHOT 
WRITE( 9 , * ) ' ANUCLD =',ANUCLD 

IF( ITER. EQ. LAST) WRITE( 9 , * ) ' RA= ' , RA , ' PR=' , PR, ' ZL= ' , ZL 
WRITE(9, * ) ' **>r>cAAA*A*****************yicyic****A***yic**' 

ENDIF 

IF( ITER. NE. LAST) RETURN 

C GET A FIELD PRINTOUT AFTER ITERATIONS STOP 
CALL PRINT 

C WRITE IN ORDER TO COMMENCE NEW PROBLEM 
REWIND( 4 ) 

WRITE ( 4) X,Y,Z,XU,YV,ZW,U,V,W,P,T 
RETURN 

QAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

ENTRY DIFFUS 
IF (NF.EQ.4) RETURN 
DO 500 1=1, LI 
DO 500 J=l,Ml 
DO 500 K=1,N1 

C DIFFUSIVITY FOR THE U, V OR W EQUATIONS 
GAM( I , J , K) =PR 

C SPECIFY ZERO DIFFUSIVITY FOR THE ZERO GRADIENT CONDITION AT Z=ZL 
C FOR THE U AND V VELOCITIES. 

IF( NF. NE. 3 )GAM( I , J, Nl )=0 . 0 
C DIFFUSIVITY FOR THE ENERGY EQUATION 
IF (NF.EQ.5) THEN 
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GAM( I , J , K ) =1 . 0 

SPECIFY ZERO DI FFUS I VI T I ES FOR THE ADIABATIC BOUNDARIES 
GAM( I , 1 , K) =0 . 0 
GAM( I ,M1 , K)=0 .0 
GAM( I , J , Nl )=0 . 0 
GAM( I, J,1)=0.0 
ENDI F 

500 CONTINUE 

IF (NF.NE.2) RETURN 

SOURCE TERMS ARE EVALUATED ONLY FOR INTERIOR POINTS 
DO 501 1=2, L2 
DO 501 J=3,M2 
DO 501 K=2,n2 

INTERPOLATE TO GET THE VALUE OF TEMPERATURE AT THE CONTROL VOLUME 
INTERFACE, SINCE THE VELOCITY U IS EVALUATED AT THE INTERFACE. 
TM=FY( J ) *T( I , J , K ) +FYM( J ) *T( I , J-1 , K ) 

CONl=RA*PR*TM 

SUPPLY ONLY PART OF THE SOURCE TERM IN THE MOMENTUM EQ . TO AVOID 
DIVERGENCE BEFORE 200 ITERATIONS 

CON( I , J , K )=FLOAT( ITER**2 ) * CON 1/2 0 0 . * *2 
IF ( ITER.GT . 200 )CON( I , J , K ) =CONl 

501 CONTINUE 

RETURN 

END 
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APPENDIX C - Natural convection in a rectangular lx>x. 



46 



The problem being considered is that of natural convection in a fluid- 
filled rectangular box. Two vertical opposite walls are maintained at a hot 
and a cold temperature as shown in Fig. C-1. The remaining box walls are 
perfectly insulated. The objective is to numerically solve for the fiuid flow 
and heat transfer inside the box arising due to buoyancy effects. The 
governing equations for the steady-state problem assuming laminar flow, 
constant properties for the fluid, no viscous dissipation and the Boussinesq 
approximation to be true are; 



Continuity 



0U , 3v ^ 3w _ 

3y 



(C-l) 



Momentum (x direction) 

D(uu) d(vu) d(wu) _ 
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Adiabatic Wall (z=0) 




Cold Wall 
(x=1) 



Fig. C-1 . Natural convection in a rectangular enclosure. 
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Momentum (z direction) 



0(u w) 
dx 



+ 



djvWj 

dy 



+ 



d(w W) 

dZ 



dZ 



+ Pr(^ + 

dx^ 



dy^ 



^2 

0 W, 






(C-4) 



Energy 

mi . mi , _^_01 ^ 

0x Dy Dz 

3.^0 + d^l , 0^0 

Dx^ dy^ dz^ (C-5) 

The above non-dimensional equations have been derived using the scaling 
of H and o/H for the lenght and velocity, respectively; where H is the 
characteristic dimension of the box in the x direction and a is the thermal 
diffusivity of the fluid. Pr=v/a is the Prandtl number of the fluid and 
Ra=gP(Th-Tc)H4/av is the Rayleigh number; v is the kinematic viscosity, g is 
the gravitational acceleration, P is the coefficient of thermal expansion and 
Th and T,, are the temperatures of the hot and cold walls, respectively. The 
temperature has been non-dimensionalized as 0=(T-Tc)/(T}^-'rj.). 

As mentioned earlier, the equations to be solved are cast into the 
canonical form in Eq. (1). The appropriate definitions for p, p, P, Sp and Sc 

are given in Table C-1: 
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Table C-1 Definition of variables in the canonical equation (Eq. 1)D) 



Equation(2) 




r 


Sp 


S (3) 


X momentum 


u 


Pr 


0 


0 


y momentum 


V 


Pr 


0 


PrRaG 


z momentum 


w 


Pr 


0 


0 


Energy 


0 


1 


0 


0 



(1) Value of p to be used in the canonical form is 1, and is set in variable 
RHOCON in the BLOCK DATA. 

(2) It is not required to cast the continuity equation in the canonical form. 

(3) The pressure gradient terms are incorporated internally in the 
program code and need not be specified as source terms in the subroutine 
USE. 



Taking advantage of the symmetry plane z=l, the computational domain is 
restricted to z- 1. The appropriate boundary conditions then become: 



ae 

U V w — 

U, V, w, 



= 0 



at z=0 



3u dv 



dd 



dz’ dz’^’dz-^ 



at z=l 



u=v=w=0, 0=1 atx=0 

U,v,w,0=O atx=l 
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U,V,W,;^= 0 



at y=0 and 1 



Program Implementation: 

As mentioned earlier, the user needs set up only the subroutine USE and 
BLOCK DATA for a particular problem to be solved. The Subroutine USE 
and other BLOCK DATA options to be used for this problem arc given at the 
end of the source code in Appendix B. The subroutine is commented 
appropriately for the user to understand the program usage. Special 
considerations in the subroutine USE are as follows: 

(1) Whenever a gradient of a dependent variable is specified to bo zero at 
a boundary, the boundary value off associated with that is set to zero in 
ENTRY DIFFUS. Thus the dependency of the boundary value of the 
variable on the solution of the variable at the adjacent node is dropped, 
leading to much faster convergence. The same solution is obtained if the 
boundary F values were kept as those in Table C-1, however, convergence 
will be slower. 

(2) For the y momentum equation, the full source term RaPrG is not 
used right from the start of the iterations, but only a portion of it is used and 
the source term is gradually increased with iterations to its fullest value. 
This helped in obtaining convergence, which could not be obtained earlier 
with the full source term for all iterations. The user is referred to Patankar 
[1], Ch. 7, for a detailed discussion on convergence and underrelaxation 
techniques for the source terms. 

(3) Since the u velocity U(I,J,K) is evaluated on the staggered grid at the 
grid-point location of the interface of the control volumes around 
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temperature 0 must be evaluated at the interface via interpolation. 

Results : 

Shown in Fig. C-2 is the local Nusselt number at the hot wall 
compared with the solution from Mallinson and Davis [2], The local 
Nusselt number has been defined as: 



As can be seen, the agreement is quite good. The small disagreement is 



discretization procedure and difference in numerical evaluation of the 
Nusselt number. 




(C.6) 



attributed to the difference in the numerical procedure, difference in 
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Local Nusselt Number 



6 





0.0 0.2 0.4 0.6 0.8 1.0 

z/ZL 



Fig. C-2. Comparison of local Nusselt numbers. 



53 



r^:feri^:nces 



1. Patankar, S. V., 1980, Numerical Heat Transfer and Fluid Flow, 
Hemisphere. 

2. Mallinson, G. D. and Davis, G. De Vahl, 1977, "Three-Dimensional 
Natural Convection in a Box: A Numerical Study," Journal of Fluid 
Mechanics, Vol. 83, Part 1, pp. 1-31. 



54 



DISTRIBUTION LIST 



1. Naval Weapons Support Center 5 

Code 6042 

Attn: Mr. Tony Buechler 
Crane, IN 47522 

2. Defense Technical Information Center 2 

Cameron Station 

Alexandria, VA 22304-6145 

3. Library, Code 0142 2 

Naval Postgraduate School 

Monterey, CA 93943-5002 

4. Superintendent 10 

Naval Postgraduate School 

Attn: Professor Y. Joshi, Code ME/Ji 
Department of Mechanical Engineering 
Monterey, CA 93943-5004 

5. Superintendent 1 

Naval Postgraduate School 

Attn: Professor M. D. Kelleher, Code ME/Kk 
Department of Mechanical Engineering 
Monterey, CA 93943-5004 

6. Research Office 1 

Naval Postgraduate School 
Code 08 

Monterey, CA 93943 



I 



I 




